Chance-Constrained DC Optimal Power Flow Using Constraint-Informed Statistical Estimation

TL;DR

This paper introduces a new method for chance-constrained DC optimal power flow that reduces complexity by focusing on aggregate and line errors, improving accuracy and efficiency over traditional high-dimensional distribution modeling.

Contribution

It proposes a constraint-informed uncertainty estimation approach that simplifies the modeling process by targeting lower-dimensional error distributions, enhancing optimization performance.

Findings

Significant reduction in modeling complexity.

Improved statistical accuracy over existing methods.

Enhanced optimization results on real-world data.

Abstract

Chance-constrained optimization has emerged as a promising framework for managing uncertainties in power systems. This work advances its application to the DC Optimal Power Flow (DC-OPF) model, developing a novel approach to uncertainty modeling and estimation. Current methods typically tackle these problems by first modeling random nodal injections using high-dimensional statistical distributions that scale with the number of buses, followed by deriving deterministic reformulations of the probabilistic constraints. We propose an alternative methodology that exploits the constraint structure to inform the uncertainties to be estimated, enabling significant dimensionality reduction. Rather than learning joint distributions of net-load forecast errors across units, we instead directly model the one-dimensional aggregate system forecast error and two-dimensional line errors weighted by power transfer distribution factors. We evaluate our approach under both Gaussian and non-Gaussian distributions on synthetic and real-world datasets, demonstrating significant improvements in statistical accuracy and optimization performance compared to existing methods.

Tables

Table 1. TABLE I: Fitted estimates of Ω \Omega on the dataset in Figure  2

Classical	${\mathrm{𝟏}}^{\top }{\widehat{𝝁}}_k$	${\mathrm{𝟏}}^{\top }{\widehat{𝚺}}_k\mathrm{𝟏}$	${\hat{\pi }}_k$
Component 1	0.24	1374.00	0.9994
Component 2	5366.81	$\approx 0.00$	0.0002
Component 3	-1760.55	407031.26	0.0004
Constraint-Informed	${\hat{m}}_k$	${\hat{\sigma }}_k^2$	${\hat{\beta }}_k$
Component 1	-0.01	4.22	0.9079
Component 2	5366.81	$\approx 0.00$	0.0002
Component 3	-4.89	30149.67	0.0919
Gaussian MLE	0.61	8536.01	1.0000

Table 2. TABLE II: Classical versus Constraint-informed tradeoffs

Approach	Number of Model Fittings	Model Dimension
Classical	$1$	$\left|ℬ\right|$
Constraint-informed	$1+\left|ℒ\right|$	$\le 2$

Table 3. TABLE III: Best BIC scores and optimal number of components K K for the aggregate forecast error Ω \Omega (lower BIC is better)

Dataset	Classical	Constraint-Informed
Synthetic-G	$-14761.09\left(K=1\right)$	$-14716.09\left(K=1\right)$
Synthetic-C	$44009.11\left(K=9\right)$	$8986.49\left(K=4\right)$
NordPool	$5721.74\left(K=9\right)$	$5517.40\left(K=6\right)$

Table 4. TABLE IV: Computational time comparisons on the IEEE 118-Bus System

	Classical	Constraint-Informed
Statistical Fitting	2.81s for $𝝃$	2.61s for $\Omega$ and 0.09s per ${𝜼}_l$
Optimization	0.34s	0.43s

Table 5. TABLE V: Computational time comparisons on the Polish 2746-bus system

	Classical	Constraint-Informed
Statistical Fitting	12.27s for $𝝃$	3.01s for $\Omega$ and 0.13s per ${𝜼}_l$
Optimization	183.97s	101.68s

Equations

$$\Omega :=1^{\top }\mathbf{\xi }=\sum _{i\in \mathcal{B}}{\xi }_i.$$

$$p_g(\mathbf{\xi }):={\bar{p}}_g-{\alpha }_g\Omega .$$

$$p_i(\mathbf{\xi })=\sum _{g\in {\mathcal{G}}_i}({\bar{p}}_g-{\alpha }_g\Omega )+({\bar{p}}_i^{\text{wind}}+{\xi }_i)-{\bar{d}}_i.$$

$$p_i^0:=\sum _{g\in {\mathcal{G}}_i}{\bar{p}}_g+{\bar{p}}_i^{\text{wind}}-{\bar{d}}_i.$$

$$\Delta p_i(\mathbf{\xi }):=-\sum _{g\in {\mathcal{G}}_i}{\alpha }_g\Omega +{\xi }_i.$$

$$p_i(\mathbf{\xi })=p_i^0+\Delta p_i(\mathbf{\xi }).$$

$$f_l(\mathbf{\xi })=\sum _{i\in \mathcal{B}}{H}_{li}{p}_i(\mathbf{\xi }).$$

$$f_l^0:=\sum _{i\in \mathcal{B}}{H}_{li}{p}_i^0.$$

$$\Delta f_l(\mathbf{\xi }):=\sum _{i\in \mathcal{B}}{H}_{li}\Delta p_i(\mathbf{\xi }).$$

$$\Delta f_l(\mathbf{\xi })=(-\sum _{i\in \mathcal{B}}{H}_{li}\sum _{g\in {\mathcal{G}}_i}{\alpha }_g)\Omega +\sum _{i\in \mathcal{B}}{H}_{li}{\xi }_i.$$

$${\gamma }_l(\mathbf{\alpha }):=-\sum _{i\in \mathcal{B}}{H}_{li}\sum _{g\in {\mathcal{G}}_i}{\alpha }_g=-({\mathbf{h}}_l^{\text{gen}}{)}^{\top }\mathbf{\alpha }.$$

$${\Lambda }_l:=\sum _{i\in \mathcal{B}}{H}_{li}{\xi }_i=({\mathbf{h}}_l^{\text{wind}}{)}^{\top }\mathbf{\xi }.$$

$$f_l(\mathbf{\xi })=f_l^0+\Delta f_l(\mathbf{\xi })=f_l^0+\left[\begin{array}{cc}{\gamma }_l(\mathbf{\alpha })& 1\end{array}\right]\left[\begin{array}{c}\Omega \\ {\Lambda }_l\end{array}\right].$$

$${\mathbf{\eta }}_l:=(\Omega ,{\Lambda }_l)=(1^{\top }\mathbf{\xi },({\mathbf{h}}_l^{\text{wind}}{)}^{\top }\mathbf{\xi })\in {\mathbb{R}}^2.$$

$$\underset{\bar{\mathbf{p}},\mathbf{\alpha }}{\min }$$

$$\sum _{g\in \mathcal{G}}{\alpha }_g=1,\mathbf{\alpha }\ge 0,\bar{\mathbf{p}}\ge 0,$$

$$\sum _{g\in \mathcal{G}}{\bar{p}}_g+\sum _{i\in \mathcal{B}}{\bar{p}}_i^{\text{wind}}-\sum _{i\in \mathcal{B}}{\bar{d}}_i=0,$$

$$\mathbb{P}\left(p_g(\mathbf{\xi })\ge p_g^{\min }\right)\ge 1-ϵ,$$

$$\mathbb{P}\left(p_g(\mathbf{\xi })\le p_g^{\max }\right)\ge 1-ϵ,$$

$$\mathbb{P}\left(f_l(\mathbf{\xi })\ge -f_l^{\max }\right)\ge 1-ϵ,$$

$$\mathbb{P}\left(f_l(\mathbf{\xi })\le f_l^{\max }\right)\ge 1-ϵ,$$

$$\sum _{g\in \mathcal{G}}{c}_{2,g}{\left({\bar{p}}_g-{\alpha }_g\mathbb{E}[\Omega ]\right)}^2+c_{2,g}{\alpha }_g^2\mathbb{V}[\Omega ]+c_{1,g}({\bar{p}}_g-{\alpha }_g\mathbb{E}[\Omega ]).$$

$$\mathbb{P}(\mathbf{a}(\mathbf{z}{)}^{\top }\mathbf{A}\mathbf{\xi }\le b(\mathbf{z}))\ge 1-ϵ,$$

$$\{{\mathbf{\xi }}^{(n)}{\}}_{n=1}^N,{\mathbf{\xi }}^{(n)}:=({\xi }_1^{(n)},{\xi }_2^{(n)},\cdots ,{\xi }_{|\mathcal{B}|}^{(n)})\in {\mathbb{R}}^{|\mathcal{B}|},$$

$$p_g(\mathbf{\xi })\sim \sum _{k=1}^K{\hat{\pi }}_k\mathcal{N}({\bar{p}}_g-{\alpha }_g{1}^{\top }{\hat{\mathbf{\mu }}}_k,{\alpha }_g^2{1}^{\top }{\hat{\mathbf{\Sigma }}}_{k}1),$$

$$f_l(\mathbf{\xi })\sim \sum _{k=1}^K{\hat{\pi }}_k\mathcal{N}(E_{lk}(\mathbf{\alpha }),V_{lk}(\mathbf{\alpha })),$$

$$E_{lk}(\mathbf{\alpha }):=f_l^0+({\gamma }_l(\mathbf{\alpha })1+{\mathbf{h}}_l^{\text{wind}}{)}^{\top }{\hat{\mathbf{\mu }}}_k,$$

$$V_{lk}(\mathbf{\alpha }):=({\gamma }_l(\mathbf{\alpha })1+{\mathbf{h}}_l^{\text{wind}}{)}^{\top }{\hat{\mathbf{\Sigma }}}_k({\gamma }_l(\mathbf{\alpha })1+{\mathbf{h}}_l^{\text{wind}}).$$

$$\sum _{k=1}^K{\hat{\pi }}_k\Phi \left(\frac{p_g^{\text{max}}-{\bar{p}}_g+{\alpha }_g{1}^{\top }{\hat{\mathbf{\mu }}}_k}{{\alpha }_g\sqrt{1^{\top }{\hat{\mathbf{\Sigma }}}_{k}1}}\right)\ge 1-ϵ,$$

$$\sum _{k=1}^K{\hat{\pi }}_k\Phi \left(\frac{f_l^{\max }-E_{lk}(\mathbf{\alpha })}{\sqrt{V_{lk}(\mathbf{\alpha })}}\right)\ge 1-ϵ.$$

$$\{{\Omega }^{(n)}{\}}_{n=1}^N,$$

$$\{{\mathbf{\eta }}_l^{(n)}{\}}_{n=1}^N,$$

$$\Omega \sim \sum _{k=1}^K{\hat{\beta }}_k\mathcal{N}({\hat{m}}_k,{\hat{\sigma }}_k^2),{\mathbf{\eta }}_l\sim \sum _{k=1}^K{\hat{\lambda }}_k\mathcal{N}(\widehat{{\mathbf{\nu }}_k},\widehat{{\mathbf{C}}_k}).$$

$$\sum _{k=1}^K{\hat{\beta }}_k\Phi \left(\frac{p_g^{\text{max}}-{\bar{p}}_g+{\alpha }_g{\hat{m}}_k}{{\alpha }_g{\hat{\sigma }}_k}\right)\ge 1-ϵ,$$

$$\sum _{k=1}^K{\hat{\lambda }}_k\Phi \left(\frac{f_l^{\max }-f_l^0-({\gamma }_l(\mathbf{\alpha }),1{)}^{\top }{\hat{\mathbf{\nu }}}_k}{\sqrt{({\gamma }_l(\mathbf{\alpha }),1{)}^{\top }{\hat{\mathbf{C}}}_k({\gamma }_l(\mathbf{\alpha }),1)}}\right)\ge 1-ϵ.$$

$${\bar{p}}_g-{\alpha }_g\left(1^{\top }{\hat{\mathbf{\mu }}}_1-{\Phi }^{-1}(1-ϵ)\sqrt{1^{\top }{\hat{\mathbf{\Sigma }}}_{1}1}\right)\le p_g^{\max }.$$

Full text

Nomenclature

$\mathcal{B}$

set of buses

$\mathcal{L}$

set of lines

$\mathcal{G}$

set of controllable generators

$\mathcal{D}$

set of loads (demands)

$\mathcal{G}_{i}$

set of generators ($\subseteq\mathcal{G}$) located at bus $i\in\mathcal{B}$

$\bm{H}$

$|\mathcal{L}|\times|\mathcal{B}|$ power transfer distribution factors

$\bm{h}^{\text{wind}}_{l},\bm{h}^{\text{gen}}_{l}$

$|\mathcal{B}|$- and $|\mathcal{G}|$-dim vectors of power transfer distribution factors mapping wind and controllable unit power outputs to flow along line $l$

$c_{1,g}$,$c_{2,g}$

linear/quadratic cost coefficient of generator $g$

$\bar{p}^{\text{wind}}_{i}$

total wind power forecast at bus $i$

$\xi_{i}$

wind power forecast error at bus $i$

$\bar{d}_{i}$

load forecast at bus $i$

$p_{g}^{\text{min}}$,$p_{g}^{\text{max}}$

min/max output of generator $g$

$f_{l}^{\text{max}}$

capacity of line $l$

$\alpha_{g}$

participation factor of generator $g$

$\bar{p}_{g},p_{g}(\bm{\xi})$

scheduled and actual output of generator $g$

$p_{i}^{0},\Delta p_{i}(\bm{\xi})$

nominal and deviation of injection at bus $i$

$f_{l}^{0},\Delta f_{l}(\bm{\xi})$

nominal and deviation of flow along line $l$

$p_{i}(\bm{\xi})$

actual net power injection into bus $i$

$f_{l}(\bm{\xi})$

actual flow along line $l$

$\gamma_{l}(\bm{\alpha})$

flow deviation along $l$ from generator control

$\Lambda_{l}$

flow deviation along $l$ from forecast errors

$\Omega$

random variable of total system forecast error

$\bm{\eta}_{l}$

random two-dim vector affecting line $l$ flow

$\hat{\pi}_{k},\hat{\bm{\mu}}_{k},\hat{\bm{\Sigma}}_{k}$

Gaussian mixture parameter estimates of $\bm{\xi}$

$\hat{\beta}_{k},\hat{m}_{k},\hat{\sigma}_{k}^{2}$

Gaussian mixture parameter estimates of $\Omega$

$\hat{\lambda}_{k},\hat{\bm{\nu}}_{k},\hat{\bm{C}}_{k}$

Gaussian mixture parameter estimates of $\bm{\eta}_{l}$

$\epsilon$

maximum acceptable violation probability

$\Phi$

standard normal distribution function

$\hat{\Phi},a_{s},b_{s}$

piecewise linear approximation of $\Phi$ and slope and intercept of the $s^{\text{th}}$ linear segment

I Introduction

Power system planning has become increasingly complex due to the presence of uncertain energy sources. Uncertainties arise from variability and unpredictability in the system, specifically due to forecast errors in net load and generation. These forecast errors directly affect the feasibility of deterministic planning procedures, such as those resulting from Optimal Power Flow (OPF) analysis, and must be accounted for[1].

To address these challenges, several approaches have been developed, including robust optimization [2, 3], multi-stage stochastic programming [4], and chance-constrained optimization [5, 6]. Chance-constrained optimization incorporates constraints that enforce certain decision-dependent random events, such as line overloads, to occur with probability no greater than a prespecified risk level. Within the OPF framework, this approach allows system operators to specify risk tolerances for power generation or transmission line violations.

Chance-constrained optimization problems are typically solved using one of two broad classes of methods: sample-based [7, 8] and sample-free analytical methods [9]. The former rely on generating a finite set of samples that represent possible realizations of the uncertain parameters; this allows converting the intractable probabilistic constraints into a deterministic form that can be solved by standard solvers. However, their key drawback is that a large number of samples may be required to obtain accurate solutions, making them computationally prohibitive for large-scale systems [10, 11]. On the other hand, sample-free methods reformulate the chance constraints into deterministic form using analytical approximations to the original problem, usually by assuming that the randomness follows certain probability distributions, such as Gaussian.

In fact, existing chance-constrained OPF (CC-OPF) models typically assume a Gaussian distribution for net load forecasting errors [6, 12, 13]. From a statistical perspective, the Gaussian assumption is appealing because the maximum likelihood estimation (MLE) of its parameters (i.e., mean and covariance) admits closed-form solutions, making parameter estimation computationally straightforward [14]. Additionally, from an optimization perspective, assuming Gaussian uncertainties enables tractable convex reformulations of the chance constraints as second-order cone programs (SOCPs), which can be efficiently solved with modern optimization solvers [6].

However, several statistical analyses have shown that non-Gaussian distributions, such as Weibull [15], Cauchy [16], generalized Pareto [17], (and Cauchy-like non-Gaussian distribution errors [18], [19]), may be better capable of representing wind generation power forecast errors, particularly within the shorter time horizon of OPF applications. Unfortunately, from the optimization perspective, the incorporation of non-Gaussian distributions is not as computationally straightforward as in the Gaussian case, since the resulting CC-OPF problems do not admit tractable reformulations that can be solved by existing solvers.

A promising path for modeling these non-Gaussian forecast errors is via Gaussian Mixture Models (GMMs), as they can approximate a wide class of probability distributions by mixing a sufficiently large number of Gaussian components [20]. At the same time, the analytical reformulation technique for Gaussian chance constraints has also been recently extended to the GMM setting [21, 22, 23, 24, 25, 26], enabling solution of these models using standard solvers. The key advantage of GMMs is that linear combinations of GMM-distributed random variables remain GMM-distributed, allowing chance constraints to be reformulated analytically using component-specific Gaussian parameters. In contrast, non-GMM distributions such as Weibull, Cauchy, or generalized Pareto, typically do not preserve their distributional form under linear combinations (e.g., the sum of two Weibull random variables is not Weibull), precluding direct analytical reformulation. Therefore, we use GMMs to approximate forecast errors, enabling both analytical reformulations and the ability to capture non-Gaussian behavior.

In this vein, [21] fits a GMM to the joint probability density of wind power output across multiple farms for chance-constrained economic dispatch; [23] uses a customized GMM to approximate the probability density of wind power generation for probabilistic OPF; [24] uses GMMs as input distributions in stochastic power flow and state estimation of distribution networks; [22] uses a generalized GMM to characterize wind power ramp distributions; and [25] and [26] fit GMMs to wind power generation at individual units for CC-OPF reformulation. In all these approaches, GMMs are fitted to the raw uncertainty at the individual component level, and the resulting distribution parameters are subsequently transformed to evaluate and reformulate the chance constraints.

However, a key drawback of this paradigm is that the joint distribution of component-level uncertainties is high-dimensional. In particular, the dimension of the forecasting error distribution scales with the number of buses with stochastic injections (e.g., wind farms). From a statistical standpoint, this approach presents several challenges. First, standard MLE algorithms become significantly harder for high-dimensional GMMs. Unlike in the case of Gaussians, the MLE objective function is now nonconvex, and the Expectation-Maximization (EM) algorithm typically used in this context is only heuristic and prone to converge to local solutions [27]. Second, the problem is affected by the curse of dimensionality: the number of parameters to be estimated can grow so large that historical data becomes insufficient, resulting in models that are prone to overfitting and that fail to generalize to unseen future data [28, 29]. In particular, roughly each additional $n$-dimensional Gaussian component in a GMM requires $O(n^{2})$ new parameters to be estimated.

Although dimensionality reduction techniques such as principal component analysis [30] or latent variable models can reduce this computational burden, they are problem-structure-agnostic and do not account for how the uncertainties actually propagate through the OPF constraints. As a result, these methods risk discarding spatial correlations among wind forecast errors or underestimating extreme events, potentially leading to constraint violations and unreliable operations. There is a need to develop an integrated framework that aligns statistical modeling with the structure of the optimization constraints and out-of-sample risk of solutions.

To address this gap, we propose a constraint-informed approach, where the OPF constraints are used to identify low-dimensional uncertainty structures that directly impact constraint feasibility. Our work highlights the importance of statistical fitting in CC-OPF problems. Instead of fitting a joint high-dimensional distribution of net-load forecast errors at system buses, we use the structure of the chance constraints to identify one- and two-dimensional uncertainties that directly impact constraint feasibility. This reduces statistical complexity and alleviates the difficulties associated with high-dimensional model fitting. Notably, the improvement stems not from modifying existing statistical estimation algorithms, but from identifying uncertainty representations relevant to the CC-OPF constraints. We demonstrate through synthetic and real-world forecast error data that this procedure improves estimation accuracy. We also illustrate the limitations of using the EM algorithm for statistical fitting, and how one can achieve better fits by modifying it from an OPF viewpoint.

Building on the constraint-informed framework, we also develop tractable analytical reformulations of chance constraints when the low-dimensional uncertainties follow a GMM. Through case studies, we show that it is the combination of constraint-informed dimensionality reduction and Gaussian mixture modeling that delivers significantly improved out-of-sample reliability compared to existing approaches. We find that the improvement is particularly stark when the underlying forecast errors exhibit unimodal but heavy-tailed behavior, elucidating that GMMs can also reduce out-of-sample risk in systems that may not involve multimodal data.

The rest of the paper proceeds as follows. Section II presents the model and assumptions; Section III presents the constraint-informed methodology, and provides reformulation of the chance constraints under Gaussian and GMM distributions; Section IV presents results of numerical experiments on synthetic and real-world datasets. Throughout, we use boldface letters for vectors and matrices, and normal font for scalars.

II Chance-Constrained DC-OPF Model

We consider a CC-OPF model with DC approximations that minimizes total generation costs while ensuring transmission line and generation limits are satisfied with high probability.

For simplicity, we assume that the demand vector $\bar{\bm{d}}$ is not random and that it can be estimated with high accuracy. This is because demand fluctuations often occur on a different timescale than the decision-making window for OPF (e.g., 5–15 minutes). To model uncertainty, we associate each wind unit with a continuous random variable representing its forecast error. These errors are then collected into a random vector $\bm{\xi}\in\mathbb{R}^{|\mathcal{B}|}$ so that $\xi_{i}$ represents the total wind power forecast error at bus $i\in\mathcal{B}$. The actual power output at bus $i$ is thus $\bar{p}^{\text{wind}}_{i}+\xi_{i}$, where $\bar{p}_{i}^{\text{wind}}$ is the forecast production. Note that we set $\bar{p}^{\text{wind}}_{i}=\xi_{i}=0$ if there are no wind units located at bus $i$. The (one-dimensional) aggregate system forecast error is:

[TABLE]

The system operator must determine the nominal power dispatch of controllable generators, denoted by $\bar{p}_{g}$, under perfect forecasts (i.e., $\bm{\xi}=\bm{0}$). These dispatch levels will serve as baseline generation and they must satisfy the expected power balance. When the actual wind power outputs realize in real time, the controllable generators will adjust their output to compensate for the resulting power imbalance. These adjustments are made under the Automatic Generation Control (AGC) policy [31], where each generator contributes proportionally to the aggregate system-wide forecast error $\Omega$. In line with existing literature [6], [32], we model this AGC policy by assigning a participation factor $\alpha_{g}\geq 0$ to each controllable generator $g\in\mathcal{G}$. Under this policy, the actual power output of generator $g$ is thus stochastic and equal to:

[TABLE]

When $\sum_{g\in\mathcal{G}}\alpha_{g}=1$, equation (2) ensures that the total adjustment matches system-wide imbalance. Furthermore, using (2), the net power injection into the system at bus $i\in\mathcal{B}$ is:

[TABLE]

We decompose this expression into a deterministic component and a stochastic deviation. First, we define the nominal power injection at bus $i$, under perfect forecast ($\bm{\xi}=\bm{0}$), as

[TABLE]

We can then model the stochastic deviation as

[TABLE]

This captures two effects happening at bus $i$ due to uncertainty. The first term reflects global balancing response from generators at bus $i$ to offset the aggregate system-wide forecast error $\Omega$, and the second term is simply the local wind forecast error at bus $i$. Using this notation, the net power injection under uncertainty can then be compactly expressed as

[TABLE]

The total power flow across line $l$ can then be expressed as a linear function of the net power injections using the matrix of Power Transfer Distribution Factors (PTDF) [33]:

[TABLE]

As before, we can use (4) and (5), to define the nominal flow along line $l$ under perfect forecast conditions as:

[TABLE]

The deviation of power injection $\Delta p_{i}(\bm{\xi})$ is also propagated linearly into transmission lines via the PTDF matrix. We define the resulting stochastic deviation in line flow as:

[TABLE]

This stochastic deviation measures how the line flows are affected by forecast errors, both globally and locally. Substituting (5) into (9) yields:

[TABLE]

In this expression, the coefficient of $\Omega$ represents the line flow deviation induced by global automatic control response from all controllable generators. We denote it as:

[TABLE]

The second component represents the line flow deviation induced by local wind forecast errors projected onto line $l$ via the PTDF matrix. We denote it as:

[TABLE]

These terms allow us to express the total flow along line $l$ as a linear function of global and local uncertainty components:

[TABLE]

Note that $\gamma_{l}(\bm{\alpha})$ is a deterministic decision-dependent constant. Therefore, the only source of randomness in (13) is:

[TABLE]

The one- and two-dimensional random variables, $\Omega$ and $\bm{\eta}_{l}$, will be the focus in the methodology section.

We note that the use of the PTDF matrix implicitly ensures power balance at the bus level. Therefore, our model follows existing literature (e.g., see [26, 32]) and avoids the need to explicitly include phase angles as decision variables. Indeed, the resulting model only has the nominal power generation $\bar{\bm{p}}$ and participation factors $\bm{\alpha}$ as decision variables. Once their optimal values are obtained, the voltage phase angles (which are also random variables) can be readily recovered under a specific realization of $\bm{\xi}$ by first constructing $\bm{p}(\bm{\xi})$ using (3) and then using the system admittance matrix.

We can now formulate the CC-OPF problem as follows:

[TABLE]

The objective function (15a) minimizes the expected total generation costs. It can be readily shown that it simplifies to:

[TABLE]

Note that if we assume $\xi_{i}$ are i.i.d. zero-mean Gaussian random variables with known variance, then (16) further simplifies and reduces to the objective function used in prior work [6, 12]. Constraints (15b) and (15c) together ensure power balance within the network under the affine control policy. Constraints (15d) and (15e) are chance constraints enforcing that the power generated lies within the minimum and maximum generation limits with (sufficiently high) probability $1-\epsilon$. Similarly, constraints (15f) and (15g) are chance constraints enforcing power flows along lines to remain less than the maximum line capacity with probability $1-\epsilon$. Here, $\epsilon$ is a risk parameter to be determined by the system operator. Lower values of $\epsilon$ reflect higher risk aversion, since they enforce a higher probability of constraint satisfaction. In practice, different $\epsilon$ values may be assigned to different classes of constraints to reflect their operational importance. We refer the reader to related works [34, 12, 35, 13] for further discussion on setting $\epsilon$ values; in our experiments, we use a common value across all constraints so that the comparison focuses on the estimation methodology.

III Constraint-Informed Reformulation

Let $\bm{z}=(\bar{\bm{p}},\bm{\alpha})$ denote the decision variables of (15). Each chance constraint in (15d)–(15g) can be written in the form

[TABLE]

where $\bm{a}(\bm{z})$ and $b(\bm{z})$ are affine functions of $\bm{z}$ and $\bm{A}$ is a constant matrix. For the power generation constraints (15e)–(15d), this follows from (1) and (2) with $\bm{A}=\bm{1}^{\top}\in\mathbb{R}^{1\times|\mathcal{B}|}$, so that $\bm{A}\bm{\xi}=\Omega\in\mathbb{R}$. For the line flow constraints (15g)–(15f), this follows from (8), (13), and (14) with $\bm{A}=\begin{bmatrix}\bm{1}&\bm{h}_{l}^{\mathrm{wind}}\end{bmatrix}^{\top}\in\mathbb{R}^{2\times|\mathcal{B}|},$ so that $\bm{A}\bm{\xi}=\bm{\eta}_{l}\in\mathbb{R}^{2}$.

The key observation is that the constraint (17) depends on uncertainty only through the projection $\bm{A}\bm{\xi}$, which is a low-dimensional linear function of the raw uncertainty $\bm{\xi}$. To evaluate the probability, one only needs the distribution of $\bm{A}\bm{\xi}$ and not of $\bm{\xi}$ itself. Since $\bm{A}\bm{\xi}$ is one-dimensional for power generation and two-dimensional for line flow constraints, learning its distribution is significantly more efficient and accurate than learning the full joint distribution of $\bm{\xi}\in\mathbb{R}^{|\mathcal{B}|}$.

We refer to this as the constraint-informed approach because the projection $\bm{A}$ is informed directly from the structure of the chance constraints. The matrix $\bm{A}$ characterizes precisely how the raw uncertainty $\bm{\xi}$ enters and influences each constraint. The approach applies to any linear chance-constrained model: one needs only to identify the low-dimensional projection $\bm{A}\bm{\xi}$ and fit a GMM distribution to the transformations $\{\bm{A}\bm{\xi}^{(n)}\}_{n=1}^{N}$ of the raw samples $\{\bm{\xi}^{(n)}\}_{n=1}^{N}$. Because the distribution is fit to low-dimensional data, it is more accurate and efficient to learn compared to fitting $\bm{\xi}$ directly and then transforming its fitted parameters. We accordingly refer to the constraint-informed approach as transform and then fit, in contrast to the classical fit and then transform. Fig. 1 illustrates this comparison.

For brevity, we present the detailed reformulation steps for the upper-bound constraints (15e) and (15g) only. The treatment of (15d) and (15f) is entirely analogous.

III-A Classical Approach: Fit and then Transform

The vast majority of existing approaches, which we shall collectively refer to as the classical approach, reformulate the chance constraints by first fitting a high-dimensional probability distribution to the raw forecast error data and then transforming the resulting distribution parameters to embed within the CC-OPF reformulation. To illustrate the main steps of the approach, consider $N$ observations of forecast errors:

[TABLE]

Fit a multivariate Gaussian–or more generally, Gaussian Mixture–distribution, $\sum_{k=1}^{K}\hat{\pi}_{k}\mathcal{N}(\hat{\bm{\mu}}_{k},\hat{\bm{\Sigma}}_{k})$, to the raw dataset $\{\bm{\xi}^{(n)}\}_{n=1}^{N}$. This is typically done using the Expectation-Maximization (EM) algorithm. 2. 2.

Transform the fitted parameters to obtain the distributions of the random variables, $p_{g}(\bm{\xi})$ and $f_{l}(\bm{\xi})$, which appear in chance constraints (15e)–(15f). Since these are affine functions of $\bm{\xi}$ (see (2)–(3), (7)–(13)), we obtain:

[TABLE] 3. 3.

Reformulate chance constraints (15e) and (15g) using the standard normal cumulative distribution function $\Phi$:

[TABLE]

Reformulate the above as convex conic constraints using the method presented in Section III-E.

III-B Constraint-Informed Approach: Transform and then Fit

The classical approach is an indirect way of targeting randomness in OPF constraints: it attempts to fit a high-dimensional joint distribution to $\{\bm{\xi}^{(n)}\}_{n=1}^{N}$ and hopes that the linear transformations of the estimators $\hat{\pi}_{k},\hat{\bm{\mu}}_{k},\hat{\bm{\Sigma}}_{k}$ accurately represent the true mean and variance of stochastic power generation $p_{g}(\bm{\xi})$ and line flow $f_{l}(\bm{\xi})$. To address this issue, we propose our constraint-informed approach, which directly targets the randomness $\Omega,\bm{\eta}_{l}$ that are present in constraints (15e), (15g), respectively. The main steps of the approach are:

Transform the data samples using (1) and (14) to obtain data samples for the one- and two-dimensional system-wide and line $l$ forecast errors, respectively:

[TABLE] 2. 2.

Fit low-dimensional GMMs to $\{\Omega^{(n)}\}_{n=1}^{N}$, $\{\bm{\eta}_{l}^{(n)}\}_{n=1}^{N}$:

[TABLE] 3. 3.

Use the fitted parameters to obtain the distributions of the random variables, $p_{g}(\bm{\xi})$ and $f_{l}(\bm{\xi})$, and to reformulate the chance constraints (15e) and (15g), respectively:

[TABLE]

Reformulate the above as convex conic constraints using the method presented in Section III-E.

When making dispatch decisions, system operators are primarily concerned with system-wide imbalances rather than individual unit deviations, thus directly modeling $\Omega$ is also operationally meaningful. On the other hand, the construction of data samples for $\bm{\eta}_{l}$ from raw forecast error data takes into consideration both aggregate system-wide forecast errors $\Omega$ and line-specific localized forecast errors $\Lambda_{l}$. Therefore, our constraint-informed approach efficiently captures both the global and local effects of forecast uncertainty, while avoiding high-dimensional statistical estimation. .

III-C Benefits of Constraint-Informed Approach

The classical and constraint-informed approaches generally yield different chance constraint reformulations for a $K$-component Gaussian mixture model, due to differences in the estimated statistical parameters. However, there is one special case where they coincide: namely, when we fit a Gaussian distribution to the forecast errors (where $K=1$), a model frequently used in existing work (e.g., see [6, 12, 13]).

In this case, it is well known that the classical chance constraint (21) can be written as a linear constraint:

[TABLE]

Similarly, the constraint-informed chance constraint (26) using a Gaussian model for $\Omega\sim\mathcal{N}(\hat{m}_{1},\hat{\sigma}_{1}^{2})$ simplifies to:

[TABLE]

It turns out that only in this special case of Gaussian distributions that are estimated using maximum likelihood estimation (MLE), we have $\bm{1}^{\top}\hat{\bm{\mu}}_{1}=\hat{m}_{1}$ and $\bm{1}^{\top}\hat{\bm{\Sigma}}_{1}\bm{1}=\hat{\sigma}_{1}^{2}$ so that the reformulated linear constraints (28) and (29) are identical. This is because of the structure of the optimal maximum likelihood estimator of Gaussians (see Appendix A). An analogous argument holds for the line flow chance constraints as well.

In reality, however, the true distribution can deviate significantly from Gaussianity, especially when the forecast error data is multi-modal or has heavy Cauchy-like tails [16, 19]. The example in Figure 2 and Table I shows that a Gaussian fit can yield a poor approximation even when aggregating forecast errors from as few as five wind units, each modeled as a Cauchy distribution with parameters from the literature [16]111The GMM is estimated using the sklearn.mixture Python package. For each method, we perform ten runs of EM with different initializations and select the model with the lowest Bayesian Information Criterion (BIC)..

In the general GMM setting, the constraint-informed approach becomes particularly advantageous, as it retains a low-dimensional uncertainty even when the underlying distribution is no longer Gaussian. In this case, the classical and constraint-informed approaches yield different reformulated constraints due to differences in their statistical parameters. In fact, even if $\bm{\xi},\Omega$ and $\bm{\eta}_{l}$ are all estimated using GMMs with the same number of components $K$, we have $\hat{\pi}_{k}\neq\hat{\beta}_{k}$, $\bm{1}^{\top}\hat{\bm{\mu}}_{k}\neq\hat{m}_{k}$, and $\bm{1}^{\top}\hat{\bm{\Sigma}}_{k}\bm{1}\neq\hat{\sigma}_{k}^{2}$ for any component $k$. This is illustrated in Table I, which compares the parameters of a 3-component GMM using both approaches on samples drawn from Figure 2. In particular, the dominant component of the constraint-informed GMM exhibits significantly lower variance compared to that of the classical approach, also evidenced by the distinct density curves in Figure 2. This can be attributed to the heuristic and non-global nature of the EM algorithm, which is typically used to estimate the model parameters.

Another key drawback of the classical approach is that it relies on high-dimensional statistical fitting, where the number of parameters grows quadratically in the number of wind units due to covariance matrix estimation. In contrast, the constraint-informed approach fits $\Omega$ and $\bm{\eta}_{l}$ directly to a one- and two-dimensional GMM, lowering the estimation burden.

III-D Computational Considerations

In the constraint-informed approach, $\Omega$ serves as a common uncertainty term for all $|\mathcal{G}|$ constraints of the form (15d) and (15e), but each of the $|\mathcal{L}|$ constraints of the form (15f) and (15g) requires a distinct random variable $\bm{\eta}_{l}$, necessitating a higher number of model fittings compared to the classical approach. This introduces a tradeoff: the classical approach fits a single high-dimensional distribution, while the constraint-informed approach fits two-dimensional distributions, reducing the computational complexity per estimation task but increasing the number of required fits; see Table II.

The statistical fitting can be efficiently performed using standard software packages that provide implementations of the EM algorithm. Notably, it can be done entirely off-line (e.g., day-ahead) once generator commitments are fixed. The online computational cost of solving the constraint-informed CC-OPF is identical to the classical approach, since the fitting does not interfere with real-time dispatch decisions.

III-E GMM Chance Constraint Reformulation

We use the constraint-informed estimates $\hat{\beta}_{k}$, $\hat{m}_{k}$, and $\hat{\sigma}_{k}^{2}$ for $\Omega$ to rewrite constraint (15e) as (26). Unlike the Gaussian case, there is no closed-form expression for the inverse cumulative distribution function (CDF) of a GMM. As a result, the left-hand side expression in (26) is both nonlinear and nonconvex, making direct reformulation difficult for optimization.

To address this challenge, we adapt an existing method from the literature [25] that builds a piecewise linear (PWL) approximation of the standard normal CDF $\Phi$. The key idea is to discretize $\Phi$ over the nonnegative reals into $S$ linear segments, using breakpoints $t_{0}=0<t_{1}<\cdots<t_{S-1}$. Given a user-specified tolerance $\delta>0$, these breakpoints are chosen optimally by [25, Algorithm 1] to minimize the number of segments $S$ while ensuring that $\lVert\hat{\Phi}(x)-\Phi(x)\rVert\leq\delta$ for all $x\in[0,\infty)$. Thus, $S$ is not a free parameter but is determined by the choice of $\delta$; a smaller $\delta$ yields a tighter approximation at the cost of more segments. We discuss how to set $\delta$ in Section IV. For each segment defined over the interval $[t_{s-1},t_{s}]$, we approximate $\Phi$ from below with a linear function whose slope $a_{s}$ and intercept $b_{s}$ can be calculated as follows:

[TABLE]

To handle the domain $[t_{S-1},\infty)$ where $\Phi$ asymptotes, the slope and intercept of the rightmost horizontal segment are:

[TABLE]

The resulting piecewise linear approximation $\hat{\Phi}$ is concave and can be expressed compactly as:

[TABLE]

The pointwise minimum ensures that the PWL approximation $\hat{\Phi}$ is a valid under-estimator for $\Phi$, so that the solution obtained by substituting $\Phi$ with $\hat{\Phi}$ in constraint (26) will also be feasible for the original constraint (15e). In other words, the PWL approximation always yields conservative (i.e., more reliable) solutions, regardless of the number of segments $S$.

To ensure that discretizing over the nonnegative reals will suffice, we need to impose additional constraints [25, Proposition 1] that restrict the domain of $\Phi$ in (26) to $[0,\infty)$.

[TABLE]

Depending on the value of $\hat{m}_{k}$, this constraint can lead to an infeasible optimization model. This can be corrected by enforcing $\hat{m}_{k}=0$ during fitting using a modified EM algorithm [36], effectively making constraint (34) non-binding.

Substituting $\Phi$ with its PWL approximation $\hat{\Phi}$ in (26) yields:

[TABLE]

We can multiply both sides by (the nonnegative term) $\alpha_{g}$ to eliminate nonlinearity in the constraint. Similarly, we can introduce a continuous auxiliary variable $M_{gk}^{1}$ for each component $k\in[K]$ and each $g\in\mathcal{G}$ to eliminate the minimum operator; we emphasize that no binary variables are needed.

[TABLE]

The reformulation of the line flow constraint (27) follows a similar idea to that of the power generation constraint (26), where we need a constraint analogous to (34):

[TABLE]

Convexity of the reformulation, however, requires the additional assumption that the component-specific covariance matrices of $\bm{\eta}_{l}$ share the same covariance matrix222This is supported by standard statistical software; e.g., spherical or tied covariance in the sklearn.mixture Python package. In our experiments, the tied structure achieved nearly identical log-likelihoods to unrestricted covariances while using fewer parameters, and this advantage is particularly pronounced when data are scarce.: $\hat{\bm{C}}_{k}=\tau_{k}^{2}\hat{\bm{C}}_{0}$. Relaxing this assumption for $\bm{\eta}_{l}$ would require introducing binary variables to handle arbitrary component covariances [37], destroying convexity and real-time applicability of the method. Utilizing techniques similar to (35) and (36), we obtain:

[TABLE]

Here, $\delta_{l}\geq 0$ is an auxiliary decision variable satisfying the convex second-order cone constraint:

[TABLE]

A detailed derivation can be found in [25, Proposition 1]. The resulting CC-OPF problem, with reformulated power generation and line flow limits, is a second-order cone convex optimization problem that can be solved using commercial solvers. A complete formulation can be found in Appendix B.

IV Computational Experiments

We use Julia 1.5.3 with PowerModels.jl [38] and the Gurobi solver. All runs were performed on a personal MacBook Air (M2 chip with 8-core CPU and 16GB RAM). The code and data files for reproducing the results are at https://github.com/Subramanyam-Lab/Constraint_Informed_CCOPF/.

IV-A Test Systems and Datasets

We perform experiments on a modified IEEE 118 bus test case. In particular, we replace 10 out of 19 conventional generators with wind units that don’t incur any generation cost. To assess scalability, we also consider a larger 2746-bus Polish system in Section IV-E. We set the chance constraint risk threshold to $\epsilon=0.05$. The accuracy of the piecewise linear approximation of $\Phi$ is set to $\delta=0.1\epsilon$.333This choice can be explained as follows. Since $0\leq\Phi(x)-\hat{\Phi}(x)\leq\delta$ for all $x\geq 0$, enforcing the chance constraints with $\hat{\Phi}$ in place of $\Phi$ can be shown to introduce at most $\delta$ additional conservatism compared to the original constraint (since the mixture weights sum to one). The factor of $0.1$ ensures that this conservatism is at least one order of magnitude smaller than $\epsilon$. With $\epsilon=0.05$, this gives $\delta=0.005$, resulting in $S=6$ PWL segments444Reducing $\delta$ results in more PWL segments $S$ that increase computational time but also provide a tighter approximation of the Gaussian CDF and hence a decrease in cost and solution conservatism. We find that tightening the tolerance to $\delta=0.002$ ($S=10$) and $\delta=0.001$ ($S=13$) increased solution time without any meaningful change in out-of-sample risk..

We experiment with three datasets $\{\bm{\xi}^{(n)}\}_{n=1}^{N}$ of forecast errors: synthetic Gaussian-distributed errors (Synthetic-G), synthetic Cauchy-distributed errors (Synthetic-C), and real-world NordPool data [39]. The choice of the Cauchy distribution is illustrative; the proposed framework applies to any continuous distribution since GMMs are universal approximators [20]. The forecast errors themselves need not follow a GMM and the latter is only used as a statistical approximation of the empirical data. For the synthetic cases, we generate ten independent datasets of forecast errors using distributional parameters555Gaussian: $\mu=-0.024,\sigma=0.036$ and Cauchy: $x_{0}=0,\gamma=0.02$. adapted from existing literature [26, 16]. Each dataset consists of 10,000 samples, with $N=8,000$ samples used for statistical estimation and subsequent optimization, while the remaining 2,000 for evaluating out-of-sample risk.

For the NordPool case, we use real 15-minute wind power production and intraday forecast data from seven locations in the CWE region: 50Hz, AMP, AT, FR, PL, TBW, and TTG. Among these, 50Hz, FR, and TTG include both onshore and offshore wind data, while the others have only onshore data. To ensure consistency with the synthetic cases, we treat onshore and offshore wind generation as separate units, resulting in a total of ten wind units. For each wind unit, we normalize forecast errors: $(\text{actual}-\text{forecast})/\text{actual}$. We collect 15 consecutive days (February 26 - March 12, 2025) of NordPool data, yielding 1,440 total samples, with 80% used for statistical estimation and optimization and the remaining 20% for out-of-sample risk evaluation.

IV-B Goodness-of-Fit

We assess the statistical goodness-of-fit of the classical and constraint-informed approaches through both visual and numerical comparisons. For visual comparisons, we plot the following three quantities in Figure 3.

True samples for aggregate system-wide forecast errors, $\{\Omega^{(n)}\}_{n=1}^{N}$, obtained using (23). 2. 2.

Density of $\Omega$ obtained using the classical ‘fit and then transform’ approach, see Section III-A for details. 3. 3.

Density of $\Omega$ using constraint-informed approach (25).

For numerical comparisons, we compute the Bayesian Information Criterion (BIC) scores of the GMMs estimated using both approaches for $K=1,2,\dots,10$ Gaussian components. For each $K$, we estimate the GMM using 10 different EM initializations and record the lowest BIC. The BIC balances goodness-of-fit with model complexity: a lower BIC indicates a better fit relative to model size [36]. Table III reports the lowest BIC and associated optimal $K$ for each dataset and approach. The corresponding density curves are plotted in Figure 3.

In the Synthetic-G dataset, where forecast errors are normally distributed, we find that the classical and constraint-informed approaches select $K=1$ with nearly identical BIC values and overlapping density curves, validating the MLE equivalence argument from Section III-C.

The results for the Synthetic-C dataset highlight the limitations of the classical approach. We believe this can be explained as follows: since the Cauchy distribution has undefined first and second moments, the conditions for the Central Limit Theorem do not apply. As a result, summing Cauchy-distributed errors does not approximate a normal distribution. The classical method, which fits a high-dimensional distribution and then reduces it to one dimension, overemphasizes the heavy tails by estimating a high-variance Gaussian component.

The failure of the classical approach is evidenced in the middle and right-hand side plots of Figure 3: even with more Gaussian components, the classical approach achieves substantially higher BIC scores than the constraint-informed approach, and its density curves deviate significantly from the true histograms of the Synthetic-C and NordPool datasets. In contrast, the constraint-informed approach, which directly fits a distribution to one-dimensional aggregated data, provides a much better match to the true samples. This is also supported by the BIC scores in Table III. This improvement in fitting is also achieved with a significantly smaller number of mixture components. For the Synthetic-C dataset, the constraint-informed approach requires only 4 Gaussian components, whereas the classical approach requires 9. Previous work [21] reported 7–9 Gaussian components for a GMM fit of $\bm{\xi}$ in a 39-bus system; our results both validate this observation and demonstrate that the proposed approach can achieve better accuracy with fewer components.

IV-C Out-of-Sample Risk and System Reliability

As noted in Section III-E, unconstrained EM may produce component mean estimates that violate (34) and (37), rendering the optimization model infeasible. To prevent this, all optimization-related results in this and subsequent sections use the FixedEM algorithm [36], which constrains every component mean to zero for both the classical ($\bm{\xi}$) and the constraint-informed approach ($\Omega,\bm{\eta}_{l}$). For each $\bm{\eta}_{l}$ fitting, we fix $K=3$ and fit under both tied and spherical covariance structures with 10 initializations each; because the zero-mean constraint reduces the effective model complexity, BIC does not favor additional components beyond $K=3$.

To evaluate the quality of CC-OPF solutions obtained, we use 20% of the dataset as a holdout, and denote it as $\mathcal{D}$. The empirical violation of a single chance constraint (i.e., generation or line flow limit) is then defined as [12]:

[TABLE]

If $\rho_{j}>\epsilon$, then the corresponding constraint does not meet the desired violation level, whereas $\rho_{j}\leq\epsilon$ indicates that the system remains secure as far as that particular generator or transmission line is concerned. Figure 4 shows the worst-case violation, $\max_{j\in\mathcal{C}}\rho_{j}$, averaged across ten runs on each dataset, where $\mathcal{C}$ is the set of all constraints in (15d)-(15g).

In the Synthetic-G dataset, the negative BIC score in Table III indicates that the obtained Gaussian parameters are a good fit to the forecast errors. This directly translates to the constraint violations being below the targeted $\epsilon=0.05$.

In the Synthetic-C dataset, we observe that both approaches yield worst-case constraint violations that exceed $\epsilon$. The violations in the classical approach average roughly $0.2$, whereas the constraint-informed approach achieves much smaller violations, averaging slightly below $0.1$ with lower overall variance across the ten runs.

Similar patterns can be seen in the NordPool dataset, where constraint-informed estimation consistently results in lower violations and variability compared to the classical approach. The NordPool dataset exhibits both skewness and a limited number of samples, leading to greater variability across runs even under the constraint-informed approach. We also note that constraint-informed violations are higher on NordPool than on Synthetic-C, suggesting that enforcing zero-mean components may be less effective when the underlying distribution is significantly skewed.

Overall, these results demonstrate that the constraint-informed dimensionality reduction consistently improves out-of-sample reliability: by fitting GMMs directly to the low-dimensional variables that enter each chance constraint, the EM algorithm better captures heavy-tailed (Synthetic-C) and skewed (NordPool) forecast errors, translating the improved goodness-of-fit into lower and less variable constraint violations.

IV-D Computational Time

Table IV reports GMM fitting times on the NordPool dataset using the FixedEM algorithm. For $\bm{\xi}$ and $\Omega$, we sweep $K\in\{1,\dots,10\}$ with 10 random initializations per $K$ and select the best model by BIC. For each $\bm{\eta}_{l}$, we fix $K=3$ and fit under both tied and spherical covariance structures with 10 initializations each. The per-line fitting time is consequently much smaller than that of $\Omega$, since fewer candidate models are evaluated. Due to the small size of the network, the optimization time for both approaches is relatively small. For the constraint-informed approach, the total computational time is 19.0 seconds, which is suitable for practical real-time deployment. Moreover, the per-line fits for $\bm{\eta}_{l}$ are independent and can be trivially parallelized.

IV-E Scalability

To evaluate the scalability of the constraint-informed approach, we consider the 2746-bus Polish 2003–2004 winter peak case [40]. Following [6], we simulate 20% wind penetration (of total load), evenly distributed across 50 wind units co-located with the 50 largest controllable generators. We use the Synthetic-C dataset where each wind unit is assigned Cauchy-distributed forecast errors.666Cauchy parameters: $x_{0}=0,\gamma=0.01$. The scale parameter $\gamma$ is reduced from $0.02$ in the IEEE-118 case to moderate the forecast error magnitude for this larger network.

Figure 5 plots the best BIC score as a function of the number of Gaussian components $K$ for both approaches as the uncertainty dimension increases from 10 to 50. As in Section IV-B, all goodness-of-fit comparisons in Figure 5 use classical EM to isolate the effect of dimensionality on statistical accuracy. The constraint-informed approach consistently achieves lower (better) BIC values across all $K$, and its 50-dimensional curve nearly overlaps with the classical fit at 10 dimensions. Moreover, the constraint-informed BIC stabilizes around $K=3$ or $4$ even at high dimensions, indicating that a small number of components suffices to capture the distribution of $\Omega$. In contrast, the classical approach does not exhibit a clear BIC minimum within the examined range; at $K\leq 2$ in 50 dimensions, the fitted GMM assigns negligible density to tail observations, causing a collapse in the log-likelihood and a sharp increase in BIC.

Table V reports fitting and optimization times on the Polish system using the FixedEM algorithm, following the same procedure as in Section IV-D. The optimization times are averaged across all 10 runs irrespective of infeasibility status. As the uncertainty dimension grows from 10 to 50, the classical fitting time increases substantially, whereas the constraint-informed fitting time for $\Omega$ remains nearly unchanged since it always performs univariate estimation. The main computational cost of the constraint-informed approach lies in the per-line fits for $\bm{\eta}_{l}$; however, the per-line fitting time is roughly constant regardless of network size, and these fits are independent and can be fully parallelized. Even when performed sequentially, the total fitting and optimization time is 531.0 seconds, which is still within a typical OPF dispatch window. On the optimization side, the classical approach yields infeasible models due to poor high-dimensional parameter estimates, while the constraint-informed approach produces a feasible solution within minutes. These results confirm that the constraint-informed approach scales to large systems, maintaining both statistical accuracy and computational tractability.

The out-of-sample worst-case constraint violation for the constraint-informed approach is also illustrated in Figure 4. The stability of the low-dimensional fitting indicates that the worst-case violation remains largely unchanged as the network scales from 118 buses to 2746 buses. In contrast, the classical approach yields infeasible models across all ten runs, due to estimation challenges in fitting 50-dimensional GMMs.

V Conclusions

We study the interplay between statistics and optimization in the chance-constrained optimal power flow (CC-OPF) problem under uncertainties arising from wind power forecast errors. Our proposed constraint-informed approach uses the OPF constraints to isolate the true relevant uncertainties to two sources of randomness: a one-dimensional aggregate system-wide forecast error and a two-dimensional line-specific flow error. By reducing dimensionality, our approach integrates statistical estimation with chance-constrained optimization, thus providing system operators with information to make more reliable dispatch decisions.

When paired with GMM models, the constraint-informed approach can help alleviate the effects of heavy-tailed and skewed forecast errors that are present in both synthetic and real-world NordPool datasets. In particular, constraint-informed GMM increases the estimation accuracy of aggregate system-wide forecast errors, while also allowing different lines to have specific covariance structures. The increase in estimation accuracy translates to improved out-of-sample risk, without sacrificing total computational time.

A limitation of the current approach is that it fits separate low-dimensional distributions for each constraint-relevant projection rather than constructing a single globally coherent low-dimensional model. While this does not affect the enforcement of individual (or even two-sided) chance constraints, it would need to be addressed in settings that require joint chance constraints coupling multiple generators or lines. Extending the constraint-informed framework to such joint formulations, as well as to time-coupled problems such as unit commitment, is a promising direction for future work.

Appendix A Equivalence of Classical and Constraint-Informed Gaussian MLE

Consider $N$ data samples $\bm{x}^{(j)}\in\mathbb{R}^{D}$, for $j=1,\dots,N$.

In the classical approach, maximum likelihood estimation (MLE) for the Gaussian parameters results in the estimates:

[TABLE]

In the constraint-informed approach, we instead transform the data first to obtain the lower-dimensional samples:

[TABLE]

We then use MLE to find estimates for the Gaussian parameters for $\Omega$ to obtain:

[TABLE]

Alternatively, we could have transformed the classical MLE parameters using the property of linear transformation of Gaussians: if $\bm{x}\sim\mathcal{N}(\hat{\bm{\mu}},\hat{\bm{\Sigma}})$, then $\bm{A}\bm{x}\sim\mathcal{N}(\bm{A}\hat{\bm{\mu}},\,\bm{A}\hat{\bm{\Sigma}}\,\bm{A}^{\top})$ for any matrix $\bm{A}$. Let $\bm{A}=\bm{1}^{\top}$ and use (41)–(43) to get:

[TABLE]

[TABLE]

We thus find that both approaches yield identical parameter estimates whenever a Gaussian model is used for fitting data.

Appendix B Constraint-Informed CC-OPF Reformulation

[TABLE]