Sieving out Unnecessary Constraints in Scenario Optimization with an Application to Power Systems

TL;DR

This paper introduces a method to eliminate unnecessary constraints in scenario optimization, improving efficiency, with demonstrated effectiveness in power grid management involving chance constraints on voltage limits.

Contribution

The paper proposes a novel approach to remove superfluous constraints after analyzing support constraints, enhancing scenario optimization techniques.

Findings

Effective constraint reduction demonstrated in power system case

Improved computational efficiency in optimal power flow problems

Validated approach with illustrative example and real-world application

Abstract

Many optimization problems incorporate uncertainty affecting their parameters and thus their objective functions and constraints. As an example, in chance-constrained optimization the constraints need to be satisfied with a certain probability. To solve these problems, scenario optimization is a well established methodology that ensures feasibility of the solution by enforcing it to satisfy a given number of samples of the constraints. The main theoretical results in scenario optimization provide the methods to determine the necessary number of samples, or to compute the risk based on the number of so-called support constraints. In this paper, we propose a methodology to remove constraints after observing the number of support constraints and the consequent risk. Additionally, we show the effectiveness of the approach with an illustrative example and an application to power distribution grid management when solving the optimal power flow problem. In this problem, uncertainty in the loads converts the admissible voltage limits into chance-constraints.

Tables

Table 1. (a) ϵ = 0.05 , d = 30 formulae-sequence italic-ϵ 0.05 𝑑 30 \epsilon=0.05,d=30

case	exact maximum value	II-A (3) $\left(N=923\right)$	II-B (5) $(N=1535,$ $r=5)$	III (7) $(N=923,$ $r=5)$
$\delta \sim 𝒩(3,1)$	$1.35$	$-0.21$	$0.31$	$0.49$
$\delta \sim 𝒰(0,1)$	$0.05$	$0.001$	$0.004$	$0.007$

Table 2. (a) ϵ = 0.05 , d = 30 formulae-sequence italic-ϵ 0.05 𝑑 30 \epsilon=0.05,d=30

case	exact maximum value	II-A (3) $\left(N=923\right)$	II-B (5) $(N=1535,$ $r=5)$	III (7) $(N=923,$ $r=5)$
$\delta \sim 𝒩(3,1)$	$1.35$	$-0.21$	$0.31$	$0.49$
$\delta \sim 𝒰(0,1)$	$0.05$	$0.001$	$0.004$	$0.007$

Table 3. (b) ϵ = 0.05 , d = 100 formulae-sequence italic-ϵ 0.05 𝑑 100 \epsilon=0.05,d=100

case	exact maximum value	II-A (3) $\left(N=2230\right)$	II-B (5) $(N=4920,$ $r=17)$	III (7) $(N=2230,$ $r=17)$
$\delta \sim 𝒩(3,1)$	$1.35$	$-0.46$	$0.31$	$0.58$
$\delta \sim 𝒰(0,1)$	$0.05$	$0.0004$	$0.004$	$0.008$

Equations

$$CCP:x^{\ast }=\arg \underset{x\in \mathcal{X}\subseteq {\mathbb{R}}^d}{\min }{c}^{T}x\text{ s.t. }P(\delta \in \Delta :x\in {\mathcal{X}}_{\delta })\ge 1-ϵ,$$

$$S{P}_N:x_N^{\ast }=\arg \underset{x\in \mathcal{X}\subseteq {\mathbb{R}}^d}{\min }{c}^{T}x\text{ s.t. }x\in \bigcap _{i\in \mathcal{N}}{\mathcal{X}}_{{\delta }^{(i)}},$$

$$S{P}_{N,k}:x_{N,k}^{\ast }=\arg \underset{x\in \mathcal{X}\subseteq {\mathbb{R}}^d}{\min }{c}^{T}x\text{ s.t. }x\in \bigcap _{i\in \mathcal{N}\setminus \{k\}}{\mathcal{X}}_{{\delta }^{(i)}},$$

$$P(V(x_N^{\ast })>ϵ)\le \beta :=\sum _{i=0}^{d-1}\left(\genfrac{}{}{0ex}{}{N}{i}\right){ϵ}^i(1-ϵ{)}^{(N-i)},$$

$$\begin{array}[]{rl}SP_{N,R}^{\mathcal{A}}\mathrel{\mathop{\ordinarycolon}}x_{N,R}^{*}=&\arg\min_{x\in\mathcal{X}\subseteq\mathbb{R}^{d}}c^{T}x\\[5.69046pt] \text{ s.t. }&x\in\bigcap_{i\in\mathcal{N}\setminus\mathcal{A}(\mathcal{N})}\mathcal{X}_{\delta^{(i)}}.\end{array}$$

$$P(V(x_{N,R}^{\ast })>ϵ)\le \beta :=\left(\genfrac{}{}{0ex}{}{R+d-1}{R}\right)\sum _{i=0}^{R+d-1}\left(\genfrac{}{}{0ex}{}{N}{i}\right){ϵ}^i(1-ϵ{)}^{(N-i)}.$$

$$0=\frac{\beta }{N+1}\sum _{m=k}^N\left(\genfrac{}{}{0ex}{}{N}{k}\right)(1-ϵ{)}^{m-k}-(\genfrac{}{}{0ex}{}{N}{k})(1-ϵ{)}^{N-k},$$

$$P(V(x_N^{\ast })>ϵ(|S_N^{\ast }|))\le \beta .$$

$$0=\frac{\beta }{N+1}\sum _{m=k}^N\left(\genfrac{}{}{0ex}{}{N}{k}\right)(1-ϵ{)}^{m-k}-(\genfrac{}{}{0ex}{}{N+R}{R}\left)\right(\genfrac{}{}{0ex}{}{N}{k})(1-ϵ{)}^{N-k},$$

$$P(V(x_{\tilde{N},R}^{\ast })>ϵ(|S_{\tilde{N},R}^{\ast }|,R))\le \beta .$$

$$\underset{x\in {\mathbb{R}}^d}{\min }{c}^{T}x\text{mbox}s.t.P(\delta :||x|{|}_2\le \delta )\ge 1-ϵ,$$

$$\left|V\right|=\left|V_0\right|+{\text{diag( ∣V0∣ )}}^{-1}(Z_p(P_G+P_L)+Z_q(Q_G+Q_L)),$$

$${\left|V\right|}_{\text{new}}=\left|V\right|+{\text{diag( ∣V0∣ )}}^{-1}(Z_p\Delta P_G+Z_q\Delta Q_G).$$

$$\begin{array}[]{l}\max\sum_{i}\Delta P_{G,i}+\Delta Q_{G,i}\\[5.69046pt] \text{s.t. }V_{\min}\textbf{1}\leq\mathinner{\!\left\lvert V\right\rvert}_{\text{new}}(\Delta P_{G},\Delta P_{G},\mathinner{\!\left\lvert V\right\rvert})\leq V_{\max}\textbf{1},\end{array}$$

$$\begin{array}[]{l}\max\sum_{i}\Delta P_{G,i}+\Delta Q_{G,i}\\ \text{s.t. }\\ P(\mathinner{\!\left\lvert V\right\rvert}\mathrel{\mathop{\ordinarycolon}}V_{\min}\textbf{1}\leq\mathinner{\!\left\lvert V\right\rvert}_{\text{new}}(\Delta P_{G},\Delta P_{G},\mathinner{\!\left\lvert V\right\rvert})\leq V_{\max}\textbf{1})\geq 95\%.\end{array}$$

$$V_{\min }\text{1}\le {\left|V\right|}^{(i)}+{\text{diag( ∣V0∣ )}}^{-1}(Z_p\Delta P_G+Z_q\Delta Q_G)\le V_{\max }\text{1}.$$

$$\begin{array}[]{rl}V_{\min}-\min_{i}\mathinner{\!\left\lvert V\right\rvert}_{l}^{(i)}&\leq\big{(}\text{diag( $\mathinner{\!\left\lvert V_{0}\right\rvert}$ )}^{-1}(Z_{p}\Delta P_{G}+Z_{q}\Delta Q_{G})\big{)}_{l}\\ V_{\max}-\max_{i}\mathinner{\!\left\lvert V\right\rvert}_{l}^{(i)}&\geq\big{(}\text{diag( $\mathinner{\!\left\lvert V_{0}\right\rvert}$ )}^{-1}(Z_{p}\Delta P_{G}+Z_{q}\Delta Q_{G})\big{)}_{l}.\end{array}$$

$$x_I^{\ast }=\arg \underset{x\in \mathcal{X}\subseteq {\mathbb{R}}^d}{\min }{c}^{T}x\text{ s.t. }x\in \bigcap _{i\in \tilde{\mathcal{N}}\setminus I}{\mathcal{X}}_{{\delta }^{(i)}}.$$

$${\Delta }_I^{\tilde{N}}=\{{\delta }^{\tilde{N}}\in {\Delta }^{\tilde{N}}:x_I^{\ast }({\delta }^{\tilde{N}})\notin {\mathcal{X}}_{{\delta }^{(i)}}\forall i\in I\}\subseteq {\Delta }^{\tilde{N}}.$$

$$\begin{array}[]{l}\{\delta^{\tilde{N}}\in\Delta^{\tilde{N}}\mathrel{\mathop{\ordinarycolon}}V(x_{\tilde{N},R}^{*}(\delta^{\tilde{N}}))>\epsilon(|S_{\tilde{N},R}^{*}|,R)\}\\[2.84544pt] \subseteq\bigcup_{I}\{\delta^{\tilde{N}}\in\Delta_{I}^{\tilde{N}}\mathrel{\mathop{\ordinarycolon}}V(x_{I}^{*}(\delta^{\tilde{N}}))>\epsilon(|S_{I}^{*}|,R)\},\end{array}$$

$$\begin{array}[]{l}P(V(x_{\tilde{N},R}^{*})>\epsilon(|S_{\tilde{N},R}^{*}|,R))\\[5.69046pt] \leq\sum_{I}P(\delta^{\tilde{N}}\in\Delta_{I}^{\tilde{N}}\wedge V(x_{I}^{*})>\epsilon(|S_{I}^{*}|,R))\\[5.69046pt] \stackrel{{\scriptstyle a)}}{{=}}\binom{\tilde{N}}{R}P(\delta^{\tilde{N}}\in\Delta_{I}^{\tilde{N}}\wedge V(x_{I}^{*})>\epsilon(|S_{I}^{*}|,R))\\[5.69046pt] =\binom{\tilde{N}}{R}P(V(x_{I}^{*})>\epsilon(|S_{I}^{*}|,R)\wedge x_{I}^{*}\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I),\end{array}$$

$$\begin{array}[]{l}P(V(x_{I}^{*})>\epsilon(|S_{I}^{*}|,R)\wedge x_{I}^{*}\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I)\\[5.69046pt] =P(\bigcup\limits_{k=0}^{d}\{V(x_{I}^{*})>\epsilon(k,R)\wedge x_{I}^{*}\notin\mathcal{X}_{\delta^{(i)}}\forall i\in I\wedge|S_{I}^{*}|=k\})\\[5.69046pt] \stackrel{{\scriptstyle b)}}{{=}}\sum\limits_{k=0}^{d}P(V(x_{I}^{*})>\epsilon(k,R)\wedge x_{I}^{*}\notin\mathcal{X}_{\delta^{(i)}}\forall i\in I\wedge|S_{I}^{*}|=k)\\[5.69046pt] \stackrel{{\scriptstyle c)}}{{=}}\sum_{k=0}^{d}\binom{N}{k}P(V(x_{I}^{*})>\epsilon(k,R)\wedge x_{I}^{*}\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I\\ \hskip 65.44142pt\wedge\;|S_{I}^{*}|=k\wedge S_{I}^{*}=\{1,\dots,k\}),\end{array}$$

$$\begin{array}[]{rl}x_{k}^{*}=\arg\min_{x\in\mathcal{X}}c^{T}x\text{ s.t. }&x\in\bigcap_{i=1,\dots,k}\mathcal{X}_{\delta^{(i)}}.\end{array}$$

$$\begin{array}[]{l}P(V(x_{I}^{*})>\epsilon(k,R)\wedge x_{I}^{*}\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I\wedge|S_{I}^{*}|=k\\ \hskip 12.80365pt\wedge\;S_{I}^{*}=\{1,\dots,k\})\\[5.69046pt] \stackrel{{\scriptstyle d)}}{{=}}P(V(x_{k}^{*})>\epsilon(k,R)\wedge x_{k}^{*}\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I\wedge|S_{k}^{*}|=k\\ \hskip 22.76228pt\wedge\;x_{k}^{*}\in\bigcap_{i=k+1,\dots,N}\mathcal{X}_{\delta^{(i)}})\\[5.69046pt] =\int_{(\epsilon(k,R),1]}P(x_{k}^{*}\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I\wedge x_{k}^{*}\in\bigcap_{i=k+1,\dots,N}\mathcal{X}_{\delta^{(i)}}\hskip 2.84544pt\big{\rvert}\\ \hskip 65.44142ptV(x_{k}^{*})=\upsilon\wedge|S_{k}^{*}|=k)dF_{k}(\upsilon)\\[5.69046pt] \stackrel{{\scriptstyle e)}}{{=}}\int_{(\epsilon(k,R),1]}(1-\upsilon)^{N-k}\upsilon^{R}dF_{k}(\upsilon),\end{array}$$

$$\begin{array}[]{l}P(V(x_{\tilde{N},R}^{*})>\epsilon(k,R))\\[5.69046pt] \stackrel{{\scriptstyle\eqref{eq:PFsuprem1}}}{{\leq}}\binom{N+R}{R}P(V(x_{I}^{*})>\epsilon(k,R)\wedge x_{I}^{*}\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I)\\[5.69046pt] \stackrel{{\scriptstyle\eqref{eq:PFsuprem2}\eqref{eq:PFsuprem3}}}{{=}}\binom{N+R}{R}\sum_{k=0}^{d}\binom{N}{k}\int_{(\epsilon(k,R),1]}(1-\upsilon)^{N-k}\upsilon^{R}dF_{k}(\upsilon).\end{array}$$

$$\begin{array}[]{rl}1&\geq P(V(x_{I}^{*})\geq 0\wedge x_{I}^{*}\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I)\\[5.69046pt] &\stackrel{{\scriptstyle\eqref{eq:PFsuprem2}\eqref{eq:PFsuprem3}}}{{=}}\sum_{k=0}^{\min(d,m)}\binom{m}{k}\int_{0}^{1}(1-\upsilon)^{m-k}\upsilon^{R}dF_{k}(\upsilon),\end{array}$$

$$\begin{array}[]{rl}\gamma=&\sup_{F_{k}\in\mathcal{C}}\binom{N+R}{R}\sum\limits_{k=0}^{d}\binom{N}{k}\int_{(\epsilon(k,R),1]}(1-\upsilon)^{N-k}\upsilon^{R}dF_{k}(\upsilon)\\[5.69046pt] &\text{s.t. }\sum\limits_{k=0}^{\min(d,m)}\binom{m}{k}\int_{0}^{1}(1-\upsilon)^{m-k}\upsilon^{R}dF_{k}(\upsilon)\leq 1,\;\forall m\geq 0.\end{array}$$

$$\begin{array}[]{rl}\bar{\gamma}_{M}=&\inf_{\lambda_{m}\geq 0}\sum_{m=0}^{M}\lambda_{m}\\[5.69046pt] &\text{s.t. }\binom{N+R}{R}\binom{N}{k}(1-\upsilon)^{N-k}\upsilon^{R}1_{(\epsilon(k,R),1]}(\upsilon)\\[5.69046pt] &\hskip 14.22636pt\leq\sum_{m=k}^{M}\lambda_{m}\binom{m}{k}(1-\upsilon)^{m-k}\upsilon^{R}\\[5.69046pt] &\hskip 14.22636pt\forall\upsilon\in[0,1],\;\forall k=0,1,\dots,d\end{array}$$

$$\begin{array}[]{l}\binom{N+R}{R}\sum_{k=0}^{d}\binom{N}{k}\int_{(\epsilon(k,R),1]}(1-\upsilon)^{N-k}\upsilon^{R}dF_{k}(\upsilon)\\[5.69046pt] =\sum_{k=0}^{d}\int_{0}^{1}\binom{N+R}{R}\binom{N}{k}(1-\upsilon)^{N-k}\upsilon^{R}1_{(\epsilon(k,R),1]}(\upsilon)dF_{k}(\upsilon)\\[5.69046pt] \stackrel{{\scriptstyle\eqref{eq:dualgammaM}}}{{\leq}}\sum_{k=0}^{d}\int_{0}^{1}\sum_{m=k}^{M}\lambda_{m}\binom{m}{k}(1-\upsilon)^{m-k}\upsilon^{R}dF_{k}(\upsilon)\\[5.69046pt] =\sum_{m=0}^{M}\lambda_{m}\sum_{k=0}^{\min(d,m)}\binom{m}{k}\int_{0}^{1}(1-\upsilon)^{m-k}\upsilon^{R}dF_{k}(\upsilon)\\[5.69046pt] \stackrel{{\scriptstyle\eqref{eq:primalgamma}}}{{\leq}}\sum\limits_{m=0}^{M}\lambda_{m}.\end{array}$$

$$\begin{array}[]{c}\binom{N+R}{R}\binom{N}{k}(1-\upsilon)^{N-k}1_{(\epsilon(k,R),1]}(\upsilon)\\[5.69046pt] \leq\frac{\beta}{N+1}\sum_{m=k}^{N}\binom{m}{k}(1-\upsilon)^{m-k},\end{array}$$

Full text

I Introduction

Many real-world optimization problems are full of uncertainties, but this is ignored in many cases. However, if this uncertainty is not taken into account, it may lead to suboptimal solutions, or what is even worse, solutions that do not satisfy the constraints of the problem. Therefore, it is essential to represent these uncertainties with random variables and stochastic constraints.

There are mainly two methodologies to deal with stochastic constraints: robust optimization and chance-constrained optimization. Robust optimization [2] aims at ensuring feasibility for any realization of the random variables in the stochastic constraints, and thus prevents the worst-case scenarios. However, robust optimization solutions may be too conservative in terms of performance, since improbable extreme values of the random variables may shrink the feasible region too much. In the worst case, this may even lead to problems with an empty feasible space.

On the other hand, in chance-constrained optimization the constraints do not need to be satisfied for every realization, but only satisfied with a given probability guarantee, like $90\%$, $95\%$ or $99\%$. These chance-constrained optimization problems are typically solved by turning the stochastic constraints into deterministic ones. Sometimes this can be done exactly, for example, in the case of Gaussian uncertainty with linear constraints; in other cases, a conservative approach can be used to ensure feasibility [3]. However, in many cases the problem may not be tractable due to various reasons, like not knowing the probability distribution of the random variables.

The scenario approach [4, 5, 6, 7, 8] is a data-driven optimization that aims at solving these intractable chance-constrained optimization problems by using only samples of scenarios, i.e., realizations of the random variables. These scenarios are used to subtitute the chance-constraint by a deterministic constraint, one for every scenario. In [6] a method is derived to determine the necessary number of scenarios to ensure a given probability guarantee. However, this method may result in too conservative solutions, that have an expected risk far below the risk threshold of the original chance-constrained problem. Therefore, in [7] this method is improved by allowing to violate some scenarios. Support constraints are defined as the scenarios whose removal results in an improvement of the objective function. In [6] it is proven that the scenario approach produces exact bounds, when the solution of the problem has as many support constraints as decision variables. This is called a fully-supported solution. When this is not the case, a better bound on the risk can be computed using the method in [8] as a function of the number of support constraints.

The natural question arising is then the following: when a not fully-supported solution is observed, can the optimal value of the objective function be improved, while still keeping the risk below its limit? The contribution of this paper consist in proposing an extension of [8], to sieve out scenarios like in [7], when a not fully-supported solution is found. An illustrative example is presented to show the effectiveness of the approach. Moreover, a power system application is presented. In this application, our scenario approach allows to find a feasible solution, that is closer to the true optimum than the standard scenario approach [6].

In the context of power systems, multiple uncertainties may affect the modelling of the optimization problems. An example is the optimal power flow, an optimization problem used to determine the set-points of controllable elements in an electrical grid. Uncertainties in loads and generation availability need to be taken into account in order to provide solutions satisfying the constraints. In [9, 10, 11] this is achieved assuming a known probability distribution of these uncertainties. In [12, 13, 14, 15, 16] this assumption is dropped, and the scenario approach is used as an appropriate tool to solve this problem.

The rest of the paper is structured as follows: Section II describes some relevant definitions and results established in the context of scenario optimization. Section III presents the main contribution of the paper: the method to remove constraints after observing the number of support constraints, with an example to illustrate its effectiveness. Section IV presents the power system application. Finally, Section V draws some conclusions and proposes future work.

II Definitions and results in scenario approach

Here we present some relevant definitions and results of the scenario approach [5, 6, 7, 8] for chance-constrained optimization problems, which will be necessary to present our approach.

Definition 1** (chance-constrained optimization)**

A chance-constrained problem $CCP$ can be defined as

[TABLE]

where $d$ is the size of the decision variables $x$, $\delta$ is a random variable with support $\Delta$ included in a measurable $\sigma$-algebra, $\mathcal{X}_{\delta}$ denotes the set of feasible solutions, which depends on multiple uncertainties represented through the vector of random variables $\delta$, and $\epsilon$ is the risk of violating the constraints.

As in [6], we will consider throughout this paper the case where $\mathcal{X}$ and $\mathcal{X}_{\delta}$ are closed convex sets; and we will assume that the feasibility domain of problem (1) has a nonempty interior, and that a solution always exists and is unique.

The problem (1) may be hard to solve depending on the shape of the feasible space $\mathcal{X}_{\delta}$, and the probability distribution of $\delta$. This probability may even not be known. Therefore, the scenario approach provides a technique to solve the problem using only samples of $\delta$ or equivalently $\mathcal{X}_{\delta}$.

II-A Exact feasibility in the scenario approach [6]

The scenario approach is based on enforcing the constraints for a given number of samples of $\delta$:

Definition 2** (scenario approach)**

The scenario approach problem $SP_{N}$ with $N$ samples: $\delta^{(i)},\;i\in\mathcal{N}=\{1,\dots,N\}$, can be defined as

[TABLE]

where $\mathcal{X}_{\delta^{(i)}}$ is the feasible space given that the random variable takes value $\delta^{(i)}$.

The notion of violation probability and support constraints are key to understand the scenario approach. They can be defined as follows:

Definition 3** (violation probability)**

The violation probability of a solution $x$ is $V(x)=P(\delta\in\Delta\mathrel{\mathop{\ordinarycolon}}x\notin\mathcal{X}_{\delta})$.

The solution $x_{N}^{*}(\delta^{N})$ of (2) and its violation probability $V(x_{N}^{*}(\delta^{N}))$ are random variables, since they depend on the set of scenarios sampled $\delta^{N}=(\delta^{(1)},\dots,\delta^{(N)})\in\Delta^{N}$. In an abuse of notation, we denote by $P(\delta^{N}\in\Delta^{N}\mathrel{\mathop{\ordinarycolon}}V(x_{N}^{*}(\delta^{N}))>\epsilon)=P(V(x_{N}^{*})>\epsilon)$ the probability (with respect to the samples $\delta^{N}$) of the violation probability.

Definition 4** (support constraint [5])**

Consider the problem without scenario $k$:

[TABLE]

scenario $k$ is a support constraint if $c^{T}x_{N,k}^{*}<c^{T}x_{N}^{*}$.

The set of indices of the support constraints for the solution $x_{N}^{*}$ is denoted as $S_{N}^{*}$. There can be at most $d$ support constraints [4], $|S_{N}^{*}|\leq d$, where $|\cdot|$ denotes the set cardinality. If $|S_{N}^{*}|=d$ then the problem is said to be fully-supported. These definitions allow to introduce the following theorem:

Theorem 1

[6]** Given a number of samples $N$, the probability of the violation probability $V(x_{N}^{*})$ exceeding a risk value $\epsilon$ can be bounded by the term $\beta$ as

[TABLE]

where equality holds if the problem is fully-supported.

Then, fixing the risk parameter $\epsilon$ and $\beta$, this theorem allows to find the right number of scenarios $N$ to limit the probability of $V(x_{N}^{*})$ being above $\epsilon$.

II-B Scenario approach discarding constraints [7]

In [7] a method is proposed to trade feasibility for performance by allowing a number of constraints to be violated. This allows to get rid of outlier scenarios that would shrink the feasible space too much.

Definition 5** (scenario approach discarding constraint)**

The scenario approach $SP_{N,R}^{\mathcal{A}}$ with $N$ samples, where $R$ are selected to be violated using some removal algorithm $\mathcal{A}(\cdot)$ such that $\mathinner{\!\left\lvert\mathcal{A}(\mathcal{N})\right\rvert}=R$, can be defined as:

[TABLE]

where any removal algorithm $\mathcal{A}$ could be valid. The only requirement for the following theorem, is that the constraints removed by algorithm $\mathcal{A}$ are almost surely violated. This can be achieved be testing if removed constraints are effectively violated, and if not, remove other constraints.

Theorem 2

[7]** If the solution $x_{N,R}^{*}$ almost surely (with respect to the set of samples $\delta^{N}$) violates the $R$ constraints removed by $\mathcal{A}$ out of the total number of constraints $N$; then the probability of the violation probability $V(x_{N}^{*})$ exceeding a risk value $\epsilon$ can be bounded by $\beta$ as

[TABLE]

Note that for the same $N$ the right-hand term of (5) is larger that the one in (3). Therefore, to achieve the same bound on $P(V(x_{N,R}^{*})>\epsilon)$, the approach with removals requires a larger number of samples $N$. The approach in (4),(5) remedies a drawback in the basic scenario approach in (2),(3): the probability density of $V(x_{N}^{*})$ may be concentrated close to [math] and so the expected value, $\mathbb{E}[V(x_{N}^{*})]\ll\epsilon$. On the other hand, for the optimal solution $x^{*}$ of (1) we would expect to have $V(x^{*})$ close to $\epsilon$, if that helps to improve the performance. This means that the basic scenario approach may be too conservative. The approach in (4),(5) can produce less conservative solutions with $\mathbb{E}[V(x_{N,R}^{*})]$ closer to $\epsilon$ [7, Appendix A].

II-C Risk given support constraints [8]

In [8] the authors analyze the case of not fully-supported problems. Instead of having a fixed risk $\epsilon$, it is defined as function of the number of support constraints $|S_{N}^{*}|$ of the solution $x_{N}^{*}$ of (2): $\epsilon(|S_{N}^{*}|)$. Throughout the rest of the paper, we will make the following assumtion as in [8]:

Assumption 1** (Non-degeneracy)**

For every $N$ and $R$, with probability $1$ (with respect to $\delta^{N}$) the solution of (2) and (4) with all the constraints, except removed constraints in (4), coincides with the solution where only the support constraints are kept.

With Assumption 1 in place we have the following result:

Theorem 3

[8]** Given $\beta\in(0,1)$, for any $k=1,\dots,d$ indicating the number of support constraints, the polynomial equation in $\epsilon$ given by

[TABLE]

has exactly one solution $\epsilon(k)\in(0,1)$, and we have

[TABLE]

It can be observed in [8, Figure 1] and in Figure 1 that the number of samples $N$ computed for a fully-supported problem, produces too conservative solutions in problems with a low number $k$ of support constraints, because $\epsilon(k)<\epsilon$ for a low $k$. Therefore, the approach in (6) gives a better bound on the risk than (3).

III Scenario approach discarding constraint given support constraints

On the one hand, the method in Section II-C allows to determine a better bound on the risk $\epsilon$ as a function of the number of support constraints $k$: $\epsilon(k)$. However, it does not provide a method to get a less conservative solution if the number of support constraints is low. On the other hand, the method in Section II-B allows to obtain potentially less conservative solutions, but it does not take into account the number of support constraint to get the bound on the risk. Therefore, in this section we present a method to react by removing constraints, when having a not fully-supported problem, and thus an $\epsilon(k)$ lower than the $\epsilon$ intended. This method is in essence a combination of the methods in presented in Sections II-B and II-C with quantitative guarantees.

Consider the problem $SP_{\tilde{N},R}^{\mathcal{A}}$, same as in (4), but with a total number of samples $\tilde{N}=N+R$ before removing $R$ samples. Here $\tilde{\mathcal{N}}$ denotes the set of all samples; $x_{\tilde{N},R}^{*}$ is the solution of the problem; and $S_{\tilde{N},R}^{*}$ the set of indices of its support constraints. Then, again with Assumption 1, we have the following theorem, which is the main result and contribution of this paper:

Theorem 4

Given $\beta\in(0,1)$, for any $k=1,\dots,d$ indicating the number of support constraints, and $R$ the number of removed constraints, the equation in $\epsilon$

[TABLE]

has exactly one solution $\epsilon(k,R)\in(0,1)$. If the solution $x_{\tilde{N},R}^{*}$ almost surely (with respect to the sample $\delta^{\tilde{N}}$) violates $R$ constraints; then the probability of the violation probability $V(x_{\tilde{N},R}^{*})$ exceeding a risk value $\epsilon$ can be bounded by $\beta$:

[TABLE]

See Appendix A for the proof.

Note that the polynomial expression in (7) differs from the one in (6) by the factor $\binom{N+R}{R}\geq 1$. As a consequence, for an $\epsilon(k)$ satisfying (6), we will have $\epsilon(k,R)\geq\epsilon(k)$ to satisfy (7). So the risk bound increases as constraints are removed. Then, it is possible to get solutions closer to the desired risk, when observing a low number of support constraints. In Figure 1, it can be observed that the risk of a solution with a low number of support constraints and $R$ removed constraints, can be lower than the risk of fully-supported solutions satisfying all constraints.

Remark 1** (Implementation)**

For any scenario optimization problem, once $N$ is fixed using an approach like (3) or (5), the solutions $\epsilon(k,R)$ of (7) can be computed for all $k$ and $R$. These solutions can be stored in a look-up table, so that they can be quickly accessed for any number $k$ of support constraints observed and $R$ constraints removed.

III-A Illustrative example

In this section we show the effectiveness of the approach with an illustrative example. Consider the following chance-constrained optimization problem

[TABLE]

where $\delta$ is a random variable, and $\epsilon$ is the risk threshold. This example corresponds to minimizing a linear function subject to the constraint that the variables are within a sphere of random radius, see Figure 2. Therefore, when applying scenario optimization, the sphere with the smallest radius will be the single support constraint, and thus the number of support constraints will be $k=1$ almost surely.

In Table I we compare the results for two probability distributions using the approaches described in (3), (5) and (7), in Sections II-A, II-B and III respectively. In this example, since the probability distribution of $\delta$ is known we can determine the exact maximum value of $\delta_{\epsilon}$ such that $||x||_{2}\leq\delta_{\epsilon}\implies P(||x||_{2}\leq\delta)\geq 1-\epsilon$.

First, it can be observed that the method described in Section III performs significantly better than the ones from Sections II-A and II-B. However, it still relatively far away from the maximum value $\delta_{\epsilon}$. Note that for II-A, negative values indicate that the problem would be infeasible. When increasing the dimension $d$ of the problem, while keeping the number of support constraints $k$ constant, from Table a) to b), the difference between methods further increases.

IV Power System Application

In power systems, the optimal power flow consists of solving an optimization problem to determine the set-points of controllable elements, like the power injection of distributed energy sources. Additionally, the grid constraints like voltage limits need to be satisfied, but this may be difficult, because some parameters in this optimization problem may be uncertain. An example could be distribution grid management, where typically only a few measurements are available [17], and thus the actual value of bus loads of the grid may be unknown. Here we combine the models and settings in [11] and [13] to build a simulation framework and apply our scenario approach (7) to power systems.

First, we consider a linear approximation of the voltage magnitudes corresponding to the Linear Coupled power flow model [18, Section 5]

[TABLE]

where $P_{L},Q_{L}$ are the active and reactive loads respectively, which we consider to be uncertain parameters; $P_{G},Q_{G}$ are the active and reactive generation, which we optimize; $Z_{p},Z_{q}$ are the known impedance matrices of the grid; $\mathinner{\!\left\lvert V\right\rvert}$ is the vector of voltage magnitudes for all nodes; and $\mathinner{\!\left\lvert V_{0}\right\rvert}$ is the known vector of voltage magnitudes under no load and generation ($P_{L},Q_{L},P_{G},Q_{G}=0$). Since $P_{L},Q_{L}$ are unknown, so will be $\mathinner{\!\left\lvert V\right\rvert}$.

When changing the values of the current operating point $P_{G},Q_{G}$ by increments $\Delta P_{G},\Delta Q_{G}$, we induce a change in the voltage magnitudes:

[TABLE]

The advantage of the expression (10) over (9) is that the uncertain parameters $P_{L},Q_{L}$ no longer appear in (10). The single uncertain paremeter is the vector $\mathinner{\!\left\lvert V\right\rvert}$, since during the optimization process we have $\Delta P_{L},\Delta Q_{L}=0$ [11]. Therefore, we only need to generate samples of $\mathinner{\!\left\lvert V\right\rvert}$ for the scenario optimization, see later in (13). Moreover, since $\mathinner{\!\left\lvert V\right\rvert}$ appears as an additive term in (10), this simplifies the process of identifying support scenarios and scenarios to remove, as we will see later in (14).

An objective for the optimal power flow could be to maximize the injection of distributed energy sources, and thus minimize the energy required from the grid, while satisfying the voltage limits $V_{\min},V_{\max}$:

[TABLE]

where 1 is a vector of $1$, and $\mathinner{\!\left\lvert V\right\rvert}_{\text{new}}(\Delta P_{G},\Delta P_{G},\mathinner{\!\left\lvert V\right\rvert})$ denotes $\mathinner{\!\left\lvert V\right\rvert}_{\text{new}}$ as a function of $\Delta P_{G},\Delta P_{G},\mathinner{\!\left\lvert V\right\rvert}$, by enforcing constraint (10). Since $\mathinner{\!\left\lvert V\right\rvert}$ in (10) is an uncertain parameter, the constraints in (11) will be stochastic. We would like to satisfy these stochastic constraints with a certain probability guarantee, for example $95\%$:

[TABLE]

Historical data or models of consumption can be used to obtain samples of $P_{L},Q_{L}$: $P_{L}^{(i)},Q_{L}^{(i)}$ for $i\in\{1,\dots,N\}$, and compute the samples $\mathinner{\!\left\lvert V\right\rvert}^{(i)}$ using (9). Then, using the approaches in [17],[13], the few real-time measurements available can be used to update the samples $\mathinner{\!\left\lvert V\right\rvert}^{(i)}$. With these samples, we can solve the stochastic optimization problem (12) using our scenario approach (7). Note that from (10) the constraints for each scenario $i$ can be represented as

[TABLE]

For the optimization problem with constraints (13), it is not possible to know beforehand the number of support constraints. However, it is possible to identify the most limiting scenarios for each element $(\cdot)_{l}$ of these vector constraints:

[TABLE]

Using (14), it becomes simple to identify which scenarios are support constraints and which ones will be violated if removed. Hence, we can apply our scenario approach (7) using the following steps:

Identify the necessary number of scenarios $\tilde{N}$ using (3) with $\epsilon=0.05$ and $\beta=10^{-3}$. Initialize the set of scenarios $\tilde{\mathcal{N}}$ and the set of removed constraints $I=\emptyset$. 2. 2.

Solve the scenario optimization problem (4) using the set of scenarios $\tilde{\mathcal{N}}$ and $\mathcal{A}(\tilde{\mathcal{N}})=I$. 3. 3.

Observe the number of support constraints $k$. Compute the number of scenarios to remove $R$ using (7), such that $\epsilon(k,R)\leq\epsilon<\epsilon(k,R+1)$. Identify the indices $\{i_{1},\dots,i_{R}\}$ of $R$ support constraint looking at (14), and add them to $I$: $I=I\cup\{i_{1},\dots,i_{R}\}$. 4. 4.

Repeat step 2). If the new solution has the same number of support constraints $k$ as the previous solution, and all constraints in $I$ are violated, finish here. If not, go back to step 3).

We have tested this new scenario approach (7) in the benchmark distribution grid 123-bus test feeder [19], see Figure 3. As in [11], we have allocated distributed energy sources in some nodes, see Figure 3; and a real-time current measurement at the branch connecting the feeder to the main grid, as in [13]. Load samples are constructed by aggregating load data of several households [20], and generation availability is constructed using real solar irradiation and wind speed data [11]. We choose the standard voltage limits $V_{\min}=0.95,V_{\max}=1.05$ [11]. We consider four approaches:

• •

Optimum: using the true values of the voltage magnitudes as if they were known.

• •

Expectation: using the estimated expected values and disregarding the uncertainty.

• •

Standard scenario: using the standard scenario approach described in (2),(3).

• •

New scenario: using the new scenario approach in (7).

and use them in a simulation to solve the stochastic optimal power flow (12) once every 15 minutes in the interval from 10am to 8pm of a day.

As it can be observed in Figure (LABEL:fig:sceobjective), both Optimum and Expectation allow to inject around $10-15\%$ more renewable energy than the scenario approaches Standard scenario and New scenario. Despite that, in Table LABEL:tab:resultsint it can be observed how the extreme values of Expectation are continuously violating the voltage limits: the maximum and minimum voltage magnitude values for all nodes have means beyond the limits $0.95,1.05$, even with a gap larger than two standard deviations in the case of the maximum. This is a consequence of disregarding the uncertainty in the data. On the other hand, both Standard scenario and New scenario satisfy the constraints. Moreover, note in Figure LABEL:fig:sceobjective that the New scenario injects around $1-2\%$ more energy (even up to $3\%$) than the Standard scenario, since some scenarios are removed. In this case, this improvement could imply millions of profit if taking into account the huge amount of distribution grids across the world. An intuitive explanation of what happens, is that out of many possible load conditions/scenarios in an electrical grid, there will be only a few that are actually limiting the current solution, the support constraints. Therefore, if removing some scenarios, the solution will perform better, but still satisfy the probability guarantee on the risk.

V Conclusions

In this work we have presented an extension of the scenario approach that allows to sieve out unnecessary constraints, when observing that a problem is not fully supported and thus too conservative. Our methodology determines the number of constraints to remove given the number of support constraints. This allows to improve the performance of the solution, while preserving the probabilistic guarantees. We have shown the relevancy of the approach using an illustrative example; then we have tested its effectiveness in a real-world power system application for solving the optimal power flow problem under uncertainty. As future work, it would be interesting to analyze if tighter bounds can be derived.

Appendix A Proof of Theorem 4

This proof combines elements from the proofs in [7, 8]. We start by defining the problem, where samples with indices in the set $I=\{i_{1},\dots,i_{R}\}$ are removed, with $|I|=R$:

[TABLE]

Consider also the set of samples violating the $I$ constraints:

[TABLE]

Since the solution $x_{\tilde{N},R}^{*}$ almost surely violates $R$ constraints, we know that, for some $I$, $x_{\tilde{N},R}^{*}=x_{I}^{*}$, and $S_{\tilde{N},R}=S_{I}$, where $S_{I}$ is the set of support constraints of $x_{I}^{*}$. Therefore, we have

[TABLE]

and thus

[TABLE]

where the equality $a)$ is due to all possible combinations of indices of the removed constraints $I$. Since the scenario sampling process is i.i.d, the probability is the same for any set of indices $I$ with same size.

Now we can look into this last probability term: $P(V(x_{I}^{*})>\epsilon(|S_{I}^{*}|,R)\wedge x_{I}^{*}(\delta^{\tilde{N}})\notin\mathcal{X}_{\delta^{(i)}}\;\forall i\in I)$. Without loss of generality, we fix the set of removed constraints $I=\{N+1,\dots,N+R\}$, so that the remaining scenarios are $\tilde{\mathcal{N}}\setminus I=\mathcal{N}$. As in Section II-C, we consider having a number $k$ of support constraints $S_{I}^{*}$, i.e., $|S_{I}^{*}|=k$, of problem (15):

[TABLE]

where equality $b)$ is due to having disjoint sets since the $|S_{I}^{*}|$ can only be a single value; and $c)$ is due to all possible combination of indices of the support constraints $S_{I}^{*}$ out of the remaining samples $N$ after removing $R$ from $\tilde{N}$. Again, since the scenario sampling process is i.i.d, the probability is the same for any set of indices $S_{I}^{*}$ of the same size.

Let $x_{k}^{*}$ be the solution using only the first $k$ scenarios:

[TABLE]

With Assumption 1, the solution $x_{I}^{*}$ is non-degenerate, i.e., if $S_{I}^{*}=\{1,\dots,k\}$, then $x_{I}^{*}=x_{k}^{*}$ with probability $1$, and we have

[TABLE]

where $F_{k}(\epsilon)=P(V(x_{k}^{*})\leq\epsilon\wedge|S_{k}^{*}|=k)$ is the $k$-th cumulative density function, i.e., the probability that $x_{k}^{*}$ is fully-supported and has violation probability less than $\epsilon$. The equality $d)$ can be proven similarly as in [8, 5.1 Proof of Theorem 1, proof A=B]. The equality $e)$ comes from the fact that for a fixed violation probability value $\upsilon$, the probability of satisfying the $N-k$ constraints with indices $\{k+1,\dots,N\}$ is $(1-\upsilon)^{N-k}$; and the probability of violating the $R$ constraints with indices $I$ is $\upsilon^{R}$.

Putting all together we have:

[TABLE]

Now we use additional information about $F_{k}(\epsilon)$. Similarly as in [8], we can derive some constraints for any generic number of scenarios $m\geq 0$ instead of $N$:

[TABLE]

where the term $\min(d,m)$ reflects that for a number of scenarios $m\leq d$, there can only be $m$ support constraints.

So we can derive an upper bound $\gamma$ on $P(V(x_{N}^{*})>\epsilon(s_{N}^{*}))$: $P(V(x_{N}^{*})>\epsilon(s_{N}^{*}))\leq\gamma$, by optimizing over all possible $F_{k}(\cdot)\in\mathcal{C}\;\forall k$, where $\mathcal{C}$ is the positive cone of generalized distribution functions:

[TABLE]

By truncating $m\leq M$, where $M\geq d$, the number of constraints are reduced. Then we get a new problem with less constraints, whose optimal value $\gamma_{M}$ satisfies $\gamma\leq\gamma_{M}$. Now we consider the dual problem of $\gamma_{M}$[21]:

[TABLE]

where $1_{(\epsilon(k,R),1]}(\upsilon)$ is the indicator function. The term $\upsilon^{R}$ can be canceled on both sides of the constraint, since for $\upsilon=0$ the constraints hold for any $\lambda_{m}$. By weak duality we have $\gamma\leq\gamma_{M}\leq\bar{\gamma}_{M}$. This can be verified for any feasible point $F_{k}$ of (20) and $\lambda_{m}$ of (21):

[TABLE]

Then, for any feasible point $\lambda_{m}\geq 0$ of (21), we have $\gamma\leq\sum_{m=0}^{M}\lambda_{m}$. Let us consider the case $M=N$, and the candidate solution $\lambda_{m}=\frac{\beta}{N+1}\;\forall m$. To ensure that this is a feasible point of (21), we need to find the smallest $\epsilon(k,R)$ for every $k\in\{0,\dots,d\}$ such that for all $\upsilon\in[0,1]$

[TABLE]

where this equation corresponds to plugging the candidate solution into the constraint in (21). Now if (7) has exactly one solution in $\epsilon(k,R)\in(0,1)$, with those $\epsilon(k,R)$ we get

[TABLE]

The proof that (7) has exactly one solution in $\epsilon\in(0,1)$, is analogous to the one in [8, 5.3 Proof of Theorem 2]. After adding the factor $\binom{N+R}{R}$ the steps of the proof still hold.