Distributionally Robust Chance Constrained Optimal Power Flow Assuming Unimodal Distributions with Misspecified Modes
Abstract
Chance constrained optimal power flow (CCOPF) formulations have been proposed to minimize operational costs while controlling the risk arising from uncertainties like renewable generation and load consumption. To solve CCOPF, we often need access to the (true) joint probability distribution of all uncertainties, which is rarely known in practice. A solution based on a biased estimate of the distribution can result in poor reliability. To overcome this challenge, recent work has explored distributionally robust chance constraints, in which the chance constraints are satisfied over a family of distributions called the ambiguity set. Commonly, ambiguity sets are only based on moment information (e.g., mean and covariance) of the random variables; however, specifying additional characteristics of the random variables reduces conservatism and cost. Here, we consider ambiguity sets that additionally incorporate unimodality information. In practice, it is difficult to estimate the mode location from the data and so we allow it to be potentially misspecified. We formulate the problem and derive a separationbased algorithm to efficiently solve it. Finally, we evaluate the performance of the proposed approach on a modified IEEE30 bus network with wind uncertainty and compare with other distributionally robust approaches. We find that a misspecified mode significantly affects the reliability of the solution and the proposed model demonstrates a good tradeoff between cost and reliability.
I Introduction
With higher penetrations of renewable generation, uncertainties have increasing influence on power system operation and hence need to be carefully considered in scheduling problems, such as optimal power flow (OPF). To manage the risk arising from uncertainties, different stochastic OPF approaches have been studied. Among these formulations, CCOPF has been proposed to directly control the constraint violation probability below a predefined threshold [1, 2, 3, 4, 5, 6, 7]. Traditional methods to solve chance constrained programs require knowledge of the joint probability distribution of all uncertainties, which may be unavailable or inaccurate. However, biased estimate may yield poor outofsample performance. Randomized techniques such as scenario approximation [8, 9], which provides a priori guarantees on reliability, require the constraints to be satisfied over a large number of uncertainty samples. The solutions from these approaches are usually overly conservative with high costs [10, 7]. Another popular approach is to assume that the uncertainties follow a parametric distribution such as Gaussian [4, 5, 7]. The resulting CCOPF is often easier to solve but the solution may have low reliability unless the assumed probability distribution happens to be close to the true one.
As an alternative, distributionally robust chance constrained (DRCC) OPF models do not depend on a single estimate of the probability distribution [11, 12, 13, 14, 15, 16, 17, 10]. More specifically, DRCC models consider a family of distributions, called the ambiguity set, that share certain statistical characteristics and requires that the chance constraint holds with respect to all distributions within the ambiguity set [18, 19, 20, 21]. Most existing work characterizes the ambiguity set based on moment information obtained from historical data of the uncertainty (see, e.g., [10, 11, 14, 13]). For example, a commonly adopted ambiguity set consists of all distributions whose mean and covariance agree with their corresponding sample estimates [10, 11, 13]. Many uncertainty distributions (e.g., those associated with wind forecast error) are unimodal and so, recently, unimodality has been incorporated to strengthen the ambiguity set and reduce the conservatism of DRCC models [13, 16, 17]. However, as compared to the moments, the mode location is more likely to be misspecified in samplebased estimation.
In this paper, we study a DRCC model with an ambiguity set based on moment and unimodality information with a potentially misspecified mode location. To the best of our knowledge, this paper is the first work discussing misspecification of a value related to a structural property, though others have considered misspecification of moments [19, 14, 18, 21, 10] and misspecification of distributions [12, 22]. Our main theoretical result shows that the distributionally robust chance constraints can be recast as a set of secondorder conic (SOC) constraints. Furthermore, we derive an iterative algorithm to accelerate solving the reformulation. In this algorithm, we begin with a relaxed formulation, and in each iteration, we efficiently find the most violated SOC constraint, if any, or terminate with a globally optimal solution. We apply the theoretical results to a direct current (DC) OPF problem and conduct a case study using a modified IEEE 30bus system with wind power. We compare our results (operational cost, reliability, computational time, and optimal solutions) to those obtained using four alternative ambiguity sets [16, 20, 17, 10].
The remainder of this paper is organized as follows. Section II empirically verifies the (multivariate) unimodality of wind forecast errors and explores misspecification of the mode location. The proposed DRCC model and ambiguity set are introduced in Section III and the main theoretical results are presented in Section IV. Section V includes the case studies and Section VI concludes the paper.
Ii Unimodality of Wind Forecast Errors & Error in Mean and Mode Estimates
In this section, we first empirically verify the unimodality of wind forecast error distributions using 10,000 data samples from [6, 7] with statistical outliers omitted (total probability ). The samples were generated using a Markov Chain Monte Carlo mechanism [23] based on real data that includes both hourly forecast and actual wind generation in Germany. In Fig. 1, we depict the histograms of univariate and bivariate wind forecast errors with bins. Both histograms empirically justify our assumption that the probability distribution of wind forecast errors is unimodal.
Next, we empirically evaluate the errors of mean and mode estimates (i.e., the peak location in the histogram). We randomly extract 100 groups of samples, each group containing 500 data points, from the wind forecast error data pool. For each group of samples, we estimate the mean by taking sample averages and estimate the mode by identifying the center of the highest bin in the 15bin histogram. In Fig. 2, we plot all the mean and mode estimates and the differences between them. From the left subfigure, we observe that sampling errors have larger impacts on mode estimates than on mean estimates. From the right subfigure, we observe that the mode estimate can deviate from the corresponding mean estimate in all directions. This indicates the importance of considering the misspecification of mode location in DRCC models, because the modemean deviation shows the skewness of the uncertainty. As a result, if we misspecify the mode location (e.g., by modeling a rightskewed distribution as a leftskewed one, see Section IIID for an example), then we may mistakenly relax the chance constraint and get poor outofsample performance.
Iii DRCC Formulation
Iiia General Formulation
In this paper, we consider the following physical constraint under uncertainty:
(1) 
where represents an dimensional decision variable, and and represent two affine functions of . Uncertainty represents an dimensional random vector defined on probability space with Borel algebra and probability distribution . The assumption that and are affine in is a standard assumption in existing DRCC models and consistent with the DRCC DC OPF.
IiiB Distributionally Robust Formulation
In reality, it may be challenging to access the (true) joint probability distribution . Oftentimes we may only have a set of historical data and certain domain knowledge of . In this case, we can consider the following distributionally robust chance constraint:
(3) 
Instead of assuming that takes a specific form, we consider an ambiguity set consisting of plausible candidates of . Then, we require that chance constraint (2) holds with respect to all distributions in .
IiiC Ambiguity Sets
In this paper, we consider three ambiguity sets, denoted as for , that are defined by a combination of moment and unimodality information. Precisely, we consider a generalized notion of unimodality defined as follows.
Definition III.1
(Unimodality [26]) For any fixed , a probability distribution on is called unimodal with mode if is nondecreasing in for every Borel set .
From the definition, we notice that parameterizes the “degree of unimodality.” When , the definition coincides with the classical univariate unimodality with mode . When , the density function of (if exists) peaks at the mode and is nonincreasing in any directions moving away from the mode. As , the requirement of unimodality gradually relaxes and eventually vanishes. Under Definition III.1, we define the following three ambiguity sets: Ambiguity set 1: (moment information only)
(4) 
Ambiguity set 2: (moment and unimodality, fixed mode)
(5) 
Ambiguity set 3: (moment and unimodality, misspecified mode)
(6) 
where and denote all probability distributions on with and without the requirement of unimodality respectively; and denote the first and second moments of ; and denotes a function returning the true mode location of with and representing a single mode value and a connected and compact set. The compact set can be constructed using possible mode estimates calculated from samples of historical data.
Among these three ambiguity sets, we use as a benchmark. Set is a special case of , i.e., only contains a single value . In practice, since the mode estimate is influenced by sampling errors, the mode estimates from data samples are not the same single values but distribute around a certain area. The shape of this area decides the underlying structural skewness in the uncertainty distribution. Hence, we compare and to see how misspecified mode estimates affect the DRCC problem. In this paper, we do not additionally consider misspecified moments since this topic has been wellstudied [19, 14, 18, 21] and our main results can be easily extended based on these existing works.
IiiD Numerical Example
We use a simple example to illustrate the impact of an inaccurate mode estimate. We assume random variable follows distribution . is a biased estimate of due to sampling errors. Both distributions are illustrated in Fig. 3, where each has zero mean and unit variance. However, is rightskewed with mode at and is leftskewed with mode at . Suppose that we try to reformulate . Based on the given distributions, we find from the correct distribution and from the biased distribution . In this example, we observe that a misspecified mode estimate could shrink the confidence bound by almost a half and significantly decrease the reliability of the solution to the chance constraint.
Iv Main Results
Iva Assumptions and prior results
To compute the exact reformulation of distributionally robust chance constraints with various ambiguity sets, we make the following assumptions.
Assumption IV.1
For , we assume that
Similarly, for , we assume that,
Assumption IV.2
For , we assume that . Similarly, for , we assume that , .
Both assumptions are standard in the related literature [27, 28, 29, 16]. Assumption IV.1 ensures that the corresponding . Assumption IV.2 ensures that the constraint is satisfied at the mode. Furthermore, we assume and , since in practice the uncertainties will at least be univariateunimodal.
Reformulations of (3) under and are derived in previous work.
Since parameter has an infinite number of choices, the reformulation in Theorem IV.2 also involves an infinite number of SOC constraints. Here we obtain a similar result for the generalized ambiguity set .
IvB Reformulation for
We now present the reformation with , which is based on Theorem IV.2:
(9)  
(10) 
Compared to (8), (9) is more complicated with two parameters and each with an infinite number of choices. To solve an optimization problem with (9), we propose an iterative solving algorithm given in Algorithm 1.
Note that the reformulated optimization problem in Step 1 contains only SOC constraints.
IvC Step 2 of Algorithm 1
The challenge is how to efficiently perform Step 2 of Algorithm 1. In the following, we assume , otherwise (9) is satisfied with regardless of the values of and . Next, we define the following terms
Since , we have where
(11) 
From Assumption IV.1, we have and transform (9) into
(12) 
Since the left side of (12) is not jointly convex or concave in and (see a proof in Appendix A), we can not find the global maximum value for the left side by simply checking the boundary values or stationary points. Therefore, we propose the following algorithm to efficiently find the global maximum.
We notice that for given a and if , the maximum value of equals with maximizer . Next, by taking the derivative of , we observe that is a strictly decreasing function of . Hence, we can compute and that cause to reach its boundary values by solving . Since and , we have and hence . Then we know as .
To efficiently solve these two equalities, we will use a golden section search by first solving for on and then for on . To efficiently apply a golden section search on , we need to find a finite upper bound instead of . The following lemma describes the selection of the finite upper bound and the best region to conduct the golden section search.
Lemma IV.1
If , . The golden section search of can be conducted on . If , . The search can be conducted on . The proof is given in Appendix B.
Furthermore, from Assumption IV.2, we have .
Based on the threshold values and , we divide our discussion into three cases.
Case 1: If , . Then from (12), we find
Then, we transform the above constraint into the following equivalent form:
(13) 
where . The left side of (13) is concave on . Define the derivative of the left side as We observe that as and

if , is the unique solution of within the domain ;

else if , .
Case 2: If , . Then from (12), we find
The above problem is a onedimensional problem on . We transform it into the following form:
(14) 
We observe that is differentiable on . Then, we know that the extreme value of happens at the critical points (boundary points , or such that that ). In the following numerical analysis, we present efficient ways to find which maximize the left side of (14).
Condition 1: If , is monotonically decreasing on and is concave on . Then,

if , ;

else if and , is the unique solution of within the domain .

else if , .
Condition 2: If , is monotonically increasing on and is convex on . Then,

if , is decreasing within the domain. To find , we follow the same discussions as in Condition 1;

else if and , is first decreasing and then increasing. Define where within the domain , , and . Then,

If , .

If or , or the unique solution of within the domain that maximizes .

If , or that maximizes .

If , equals the unique solution of within the domain .

If or , .


else if , is convex on . or that maximizes .
Case 3: If , . Then from (12), we find
which we transform into the following equivalent form
(15) 
where . Define the derivative of the left hand side of (15) as Then is concave on and as , . Then,

if , is an increasing function and ;

else if , as , . Based on the concavity of , we find

if , ;

else if , equals the unique solution of within the domain .

To efficiently apply the golden section search, we determine an effective finite upper bound instead of . Let the effective upper bound be , we have
Lemma IV.2
Then, instead of a search on , we only need to search on .
Combining all three cases, we can find the overall worst case and given . If (12) is satisfied with these parameters, then there is no violated constraint in Step 2 of Algorithm 1. If (12) is not satisfied, we need to use the worst case and in Step 3 and the iteration continues. Depending on how we define , are different functions of .
IvD Candidates of
In this section, we demonstrate how the selection of affects the determination of , , and . Specifically, we give two examples of and show how to exactly reformulate (10) (i.e., Assumption IV.2) and how to calculate and , given . Furthermore, we show how to find the worst case from .
Rectangular Support: We assume that and hence we can reformulate (10) as
(16) 
Furthermore, given , we have the following relationships due to (11).
(17)  
(18) 
Based on (17) and (18), if we have the worst case , we find the worst case by solving (19) for and substituting in (20):
(19)  
(20) 
where returns a diagonal matrix whose diagonal elements equal the sign of each elements in .
V Case Study
Va Simulation Setup
We consider the DC OPF problem from [16]. We assume that the system has two wind power plants with wind forecast error . With generators and buses, the design variables are generation , up and down reserve capacities , and a distribution vector , which determines the realtime reserve provision from each generator used to balance the wind forecast error. The full problem formulation is as follows.
(26a)  
(26b)  
(26c)  
(26d)  
(26e)  
(26f)  
(26g)  
(26h)  
(26i)  
(26j) 
where , , and are cost parameters. Constraint (26b) bounds the power flow, which is calculated from the power injections defined in (26d) and the parameter matrix , by the line limits . Constraint (26c) computes the realtime reserve usage for each generator. In (26d) is the wind forecast, is the load, and , , and are matrices that map generators, wind power plants, and loads to buses; (26e) restricts generation to within its limits ; (26f) restricts by the reserve capacity; (26g), (26h) enforce power balance with and without wind forecast error; and (26i), (26j) ensure all decision variables are nonnegative.
We test our approach on a modified IEEE 30bus system with network and cost parameters from [31]. We set . We add the wind power plants to buses 22 and 5 and set MW. We use the same wind power forecast uncertainty data ( scenarios) as in Section II. We congest the system by increasing each load by and reducing the limit of the line connecting buses 1 and 2 to 30 MW. All optimization problems are solved using CVX with the Mosek solver [32, 33].
To construct the ambiguity sets, unlike in Section II, the outliers are used when estimating the statistical parameters (first moment , second moment , and the set of the mode ) and evaluating the reliability of the solution. We set , , and assume is a rectangular set.
VB Additional Ambiguity Sets
We benchmark against two additional ambiguity sets from related work.
Ambiguity set 4: (moment and unimodality with fixed mode at the mean [20])
(27) 
Ambiguity set 5: (moment and unimodality with and arbitrary mode [17])
(28) 
Set is a special case of with the mode at the mean, while is a special case of with and that is an ellipsoidal set based on and as shown in Assumption IV.1. In other words, our is more general than , , and . The reformulations of and are simpler than with a single SOC constraint
(29) 
VC Simulation Results
VC1 Estimation of
We next analyze how the data size of each sample and the number of bins within the histogram affect the estimate of the mode support. Figure 4 shows that if we change from 15 to 30 the histograms no longer show a unimodal distribution, as compared to Fig. 1. The problem is exacerbated as grows.
We next explore the impact of the size of the data pool. We first use the entire data pool to select 100 samples with different data sizes ( and ) and number of bins ( and ) and show scatter plots of the mode values in Fig. 5. As gets larger, the mode values are more condensed and hence more accurate. When and mode values appear in several disjoint regions, but this disjointness is mitigated as increases to . Based on the scatter plots, we determined the parameters of the four rectangular sets used in . The results are given in Table II.
15  Plant 1  4.44  0.10  3.45  0.17 

Plant 2  4.45  0.24  3.69  0.11  
30  Plant 1  4.36  0.19  3.02  0.93 
Plant 2  4.22  0.22  3.06  0.39 
10  Plant 1  4.77  0.58  3.52  0.09 

Plant 2  5.05  0.44  4.43  0.06  
20  Plant 1  5.82  0.06  4.36  0.09 
Plant 2  5.76  0.04  3.86  0.19 
We repeated the analysis using only a partial data pool, specifically, we randomly selected 1000 data from the full pool to comprise the partial pool. We also use different choices of and . The scatter plots are shown in Fig. 6 and parameter values for are given in Table II.
VC2 Objective Costs
We next analyze the objective costs and the optimal reserve capacities using different ambiguity sets. The results are summarized in Table III. In all case studies, since we focus on mode misspecification not moment misspecification, moments are calculated using the full or partial data pool and all ambiguity sets use the same moments.
Full pool  
M  M  M  M  M  M  
Total Cost  26160  19440  19546  18993  19547  19526  19542  19949  19818  19982  19896  19818  19982 
Generation Cost  13032  11515  11504  11506  11481  11491  11506  11522  11522  11522  11522  11514  11522 
Reserve Cost  13129  7925  8042  7488  8065  8035  8036  8427  8296  8460  8373  8304  8460 
Up Reserve (MW)  38.8  26.8  26.2  26.3  25.1  26.1  26.1  27.1  26.8  27.1  27.0  26.7  27.1 
Down Reserve (MW)  26.9  12.9  14.0  11.1  15.2  14.1  14.1  15.1  14.7  15.2  14.9  14.8  15.2 
Partial pool  
M  M  M  M  M  M  
Total Cost  21845  16759  16740  15250  15915  16735  16718  17883  17909  17876  17909  16791  17949 
Generation Cost  11713  11376  11376  11324  11295  11336  11367  11420  11420  11420  11420  11371  11420 
Reserve Cost  10132  5383  5364  3926  4620  5399  5351  6463  6489  6456  6489  5420  6529 
Up Reserve (MW)  31.1  20.1  20.1  17.5  16.1  19.4  19.1  22.0  22.0  21.9  22.0  19.2  22.1 
Down Reserve (MW)  19.6  6.8  6.7  2.1  7.0  7.6  7.7  10.3  10.4  10.4  10.4  7.9  10.6 
For ambiguity set , we perform tests with the following six fixed mode estimates.

M1: mode determined using the full (partial) data pool with histogram of () bins. This case demonstrates the performance of with an accurate mode estimate.

M2: mode determined using the full (partial) data pool with histogram of () bins. This case shows how affects the result.

M36: combinations of the largest