Distributionally Robust Chance Constrained Optimal Power Flow Assuming Unimodal Distributions with Misspecified Modes

# Distributionally Robust Chance Constrained Optimal Power Flow Assuming Unimodal Distributions with Misspecified Modes

Bowen Li,  Ruiwei Jiang,  and Johanna L. Mathieu,  This work is supported by the U.S. National Science Foundation Awards CCF-1442495 and CMMI-1662774. B. Li and J. L. Mathieu are with the Department of Electrical Engineering and Computer Science, University of Michigan at Ann Arbor, Ann Arbor, MI 48109 USA (e-mail: libowen@umich.edu; jlmath@umich.edu). R. Jiang is with the Department of Industrial and Operations Engineering, University of Michigan at Ann Arbor, Ann Arbor, MI 48109 USA (e-mail: ruiwei@umich.edu).
###### Abstract

Chance constrained optimal power flow (CC-OPF) formulations have been proposed to minimize operational costs while controlling the risk arising from uncertainties like renewable generation and load consumption. To solve CC-OPF, 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 separation-based algorithm to efficiently solve it. Finally, we evaluate the performance of the proposed approach on a modified IEEE-30 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 trade-off between cost and reliability.

Optimal power flow, chance constraint, distributionally robust optimization, misspecified mode, -unimodality

## 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, CC-OPF has been proposed to directly control the constraint violation probability below a pre-defined 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 out-of-sample 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 CC-OPF 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 sample-based 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 second-order 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 30-bus 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 15-bin 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 mode-mean deviation shows the skewness of the uncertainty. As a result, if we misspecify the mode location (e.g., by modeling a right-skewed distribution as a left-skewed one, see Section III-D for an example), then we may mistakenly relax the chance constraint and get poor out-of-sample performance.

## Iii DRCC Formulation

### Iii-a General Formulation

In this paper, we consider the following physical constraint under uncertainty:

 a(x)⊤ξ≤b(x), (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.

To manage constraint violations due to uncertainty, one natural way is to ensure that (1) is satisfied with at least a pre-defined probability threshold , which leads to the following chance constraint [24, 25]:

 Pξ(a(x)⊤ξ≤b(x))≥1−ϵ, (2)

where normally takes a large value (e.g., ).

### Iii-B 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:

 infPξ∈DξPξ(a(x)⊤ξ≤b(x))≥1−ϵ. (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 .

### Iii-C 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 non-decreasing 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 non-increasing 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)

 D1ξ:={Pξ∈Pn:EPξ[ξ]=μ, EPξ[ξξ⊤]=Σ}, (4)

Ambiguity set 2: (moment and -unimodality, fixed mode)

 D2ξ:={Pξ∈Pnα∩D1ξ: M(ξ)=mt}, (5)

Ambiguity set 3: (moment and -unimodality, misspecified mode)

 D3ξ:={Pξ∈Pnα∩D1ξ: M(ξ)∈Ξ}, (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 well-studied [19, 14, 18, 21] and our main results can be easily extended based on these existing works.

### Iii-D 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 right-skewed with mode at and is left-skewed 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

### Iv-a 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

 (α+2α)(Σ−μμ⊤)≻1α2(μ−mt)(μ−mt)⊤.

Similarly, for , we assume that,

 (α+2α)(Σ−μμ⊤)≻1α2(μ−m)(μ−m)⊤.
###### 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 univariate-unimodal.

Reformulations of (3) under and are derived in previous work.

###### Theorem IV.1

(Theorem 2.2 in [30]) With , (3) can be exactly reformulated as

 √(1−ϵϵ)a(x)⊤(Σ−μμ⊤)a(x)≤b(x)−a(x)⊤μ. (7)
###### Theorem IV.2

(Theorem 1 in [16]) With , (3) can be exactly reformulated as

 √1−ϵ−τ−αϵ∥Λta(x)∥≤τ(b(x)−μ⊤a(x)) +(τ−α+1α)(μ−mt)⊤a(x),∀τ≥(11−ϵ)1/α, (8)

where .

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 .

### Iv-B Reformulation for D3ξ

We now present the reformation with , which is based on Theorem IV.2:

 √1−ϵ−τ−αϵ∥Λa(x)∥≤(τ−α+1α)(μ−m)⊤a(x) +τ(b(x)−μ⊤a(x)), ∀τ≥(11−ϵ)1/α, ∀m∈Ξ, (9) a(x)⊤m≤b(x), ∀m∈Ξ, (10)

where and (10) comes from Assumption IV.2.

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.

### Iv-C 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

 h=a(x∗i)⊤(μ−m)/α, ~c=b(x∗i)−μ⊤a(x∗i), ~R=√a(x∗i)⊤(α+2α)(Σ−μμ⊤)a(x∗i), g(τ)=√1−ϵ−τ−αϵ, f(τ)=−(ατ−α−1).

Since , we have where

 ¯¯¯h=maxm∈Ξa(x∗i)⊤(μ−m)/α, h––=minm∈Ξa(x∗i)⊤(μ−m)/α. (11)

From Assumption IV.1, we have and transform (9) into

 [g(τ)√~R2−h2+f(τ)h]−~cτ≤0, ∀h∈[h––,¯¯¯h], ∀τ≥τ0. (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

 g(τ)√~R2−¯¯¯h2+f(τ)¯¯¯h−~cτ≤0.

Then, we transform the above constraint into the following equivalent form:

 F1(τ)=C1g(τ)−(~c+α¯¯¯h)τ+(α+1)¯¯¯h≤0, (13)

where . The left side of (13) is concave on . Define the derivative of the left side as We observe that as and

1. if , is the unique solution of within the domain ;

2. else if , .

Case 2: If , . Then from (12), we find

 ~R√g(τ)2+f(τ)2−~cτ≤0.

The above problem is a one-dimensional problem on . We transform it into the following form:

 F2(τ)=~R2(g(τ)2+f(τ)2)−~c2τ2≤0. (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).

The first and second derivative of the left side of (14) are

 F′2(τ) =~R2(αϵτ−α−1+2α(ατ−α−1))−2~c2τ =α~R2ϵτ−α−1+(2α2~R2−2~c2)τ−2~R2α(α+1), F′′2(τ) =−α~R2(α+1)ϵτ−α−2+(2α2~R2−2~c2).

Given this, there are two conditions.

Condition 1: If , is monotonically decreasing on and is concave on . Then,

1. if , ;

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

3. else if , .

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

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

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

1. If , .

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

3. If , or that maximizes .

4. If , equals the unique solution of within the domain .

5. If or , .

3. else if , is convex on . or that maximizes .

Case 3: If , . Then from (12), we find

 g(τ)√~R2−h––2+f(τ)h––−~cτ≤0,

which we transform into the following equivalent form

 F3(τ)=C3g(τ)−(~c+αh––)τ+(α+1)h––≤0, (15)

where . Define the derivative of the left hand side of (15) as Then is concave on and as , . Then,

1. if , is an increasing function and ;

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

1. if , ;

2. 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

 F′3(τ2)=C3g′(τ2)−(~c+α–h)≤0.
###### Lemma IV.2

A feasible selection of is

 τ2=[−1+√1+4(1−ϵ)C22C2]−1α,

where . The proof is given in Appendix C.

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 .

### Iv-D 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

 a(x)⊤(k––+¯¯¯k2)+|a(x)|⊤(¯¯¯k−k––2)≤b(x). (16)

Furthermore, given , we have the following relationships due to (11).

 h––=[a(x∗i)⊤(μ−–k+¯¯¯k2)−|a(x∗i)|⊤(¯¯¯k−k––2)]/α, (17) ¯¯¯h=[a(x∗i)⊤(μ−k––+¯¯¯k2)+|a(x∗i)|⊤(¯¯¯k−k––2)]/α. (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):

 h∗=[a(x∗i)⊤(μ−k––+¯¯¯k2)+λr|a(x∗i)|⊤(¯¯¯k−k––2)]/α, (19) m∗=(k––+¯¯¯k2)−λrsign(a(x∗i))(¯¯¯k−k––2), (20)

where returns a diagonal matrix whose diagonal elements equal the sign of each elements in .

Ellipsoidal Support: We assume that , where . Then we can reformulate (10) as

 a(x)⊤mc+∥∥P1/2a(x)∥∥2≤b(x). (21)

Furthermore, due to (11), we have the following relationships:

 h––=[a(x∗i)⊤(μ−mc)−∥∥P1/2a(x∗i)∥∥2]/α, (22) ¯¯¯h=[a(x∗i)⊤(μ−mc)+∥∥P1/2a(x∗i)∥∥2]/α. (23)

Next, if we have the worst case , we find the worst case directly by solving (24) for and substituting in (25):

 h∗=[a(x∗i)⊤(μ−mc)−λe∥∥P1/2a(x∗i)∥∥2]/α, (24) m∗=mc+λePa(x∗i)∥∥P1/2a(x∗i)∥∥2. (25)

## V Case Study

### V-a 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 real-time reserve provision from each generator used to balance the wind forecast error. The full problem formulation is as follows.

 min PTG[C1]PG+CT2PG+CTR(RupG+RdnG) (26a) s.t. −Pl≤APinj≤Pl (26b) RG=−dG(~w1+~w2) (26c) Pinj=CG(PG+RG)+CW(PfW+~w)−CLPL (26d) P––G≤PG+RG≤¯¯¯¯PG (26e) −RdnG≤RG≤RupG (26f) 11×NGdG=1 (26g) 11×NB(CGPG+CWPfW−CLPL)=0 (26h) PG≥0NG×1, dG≥0NG×1 (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 real-time 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 non-negative.

We test our approach on a modified IEEE 30-bus 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.

We benchmark against two additional ambiguity sets from related work.
Ambiguity set 4: (moment and unimodality with fixed mode at the mean [20])

 D4ξ:={Pξ∈Pnα∩D1ξ: M(ξ)=μ}. (27)

Ambiguity set 5: (moment and unimodality with and arbitrary mode [17])

 D5ξ:={Pξ∈Pn1∩D1ξ}. (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

 K√a(x)⊤(Σ−μμ⊤)a(x)≤b(x)−a(x)⊤μ, (29)

where can be found in [20] for and in [17] for .

### V-C Simulation Results

#### V-C1 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.

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.

#### V-C2 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.

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.

• M3-6: combinations of the largest