Robust AC Optimal Power Flow

TL;DR

This paper develops a tractable approximation method for the complex robust AC optimal power flow problem, enabling reliable operation of power systems with renewable energy uncertainties.

Contribution

It introduces a finite-dimensional semidefinite programming approach with guarantees for feasible affine recourse policies in robust AC OPF.

Findings

Effective on WSCC 9-bus system with renewable uncertainty

Applicable to IEEE 14-bus test system with high renewable penetration

Provides a sufficient condition for policy feasibility

Abstract

There is a growing need for new optimization methods to facilitate the reliable and cost-effective operation of power systems with intermittent renewable energy resources. In this paper, we formulate the robust AC optimal power flow (RAC-OPF) problem as a two-stage robust optimization problem with recourse. This problem amounts to a nonconvex infinite-dimensional optimization problem that is computationally intractable, in general. Under the assumption that there is adjustable generation or load at every bus in the power transmission network, we develop a technique to approximate RAC-OPF from within by a finite-dimensional semidefinite program by restricting the space of recourse policies to be affine in the uncertain problem data. We establish a sufficient condition under which the semidefinite program returns an affine recourse policy that is guaranteed to be feasible for the original RAC-OPF problem. We illustrate the effectiveness of the proposed optimization method on the WSCC 9-bus and IEEE 14-bus test systems with different levels of renewable resource penetration and uncertainty.

Tables

Table 1. TABLE I : (WSCC 9-bus system).   Specification of each generator’s location, marginal cost, and constraint parameters. The † † \dagger symbol indicates that the corresponding value in the table is determined by equation ( 16 ).

				Inflexible Gen.		Intermittent Gen.						Flexible Gen.			Load Buses
Parameters		Units		1		2	3	4	5	6		7	8		9
${\alpha }_i$		$/MW		30		0	0	0	0	0		50	50		–
$\left\{d_i\right\}$		MW		0		0	0	100	0	125		0	0		90
$\mathrm{Im}\left\{d_i\right\}$		MVAR		0		0	0	35	0	50		0	0		30
$\left\{g_i^{\max }\right\}$		MW		200		30	30	30	30	30		250	270		–
$\mathrm{Im}\left\{g_i^{\max }\right\}$		MVAR		300		$†$	$†$	$†$	$†$	$†$		300	300		–
$\left\{g_i^{\min }\right\}$		MW		10		0	0	0	0	0		10	10		–
$\mathrm{Im}\left\{g_i^{\min }\right\}$		MVAR		-300		$†$	$†$	$†$	$†$	$†$		-300	-300		–
$\left\{r_i^{\max }\right\}=-\left\{r_i^{\min }\right\}$		MW		0		$\infty$	$\infty$	$\infty$	$\infty$	$\infty$		$\infty$	$\infty$		–
$\mathrm{Im}\left\{r_i^{\max }\right\}=-\mathrm{Im}\left\{r_i^{\min }\right\}$		MVAR		0		$\infty$	$\infty$	$\infty$	$\infty$	$\infty$		$\infty$	$\infty$		–
Bus index (per [32])				2		4	6	7	8	9		1	3		5

Table 2. TABLE II : (IEEE 14-bus system).   Specification of each generator’s location, marginal cost, and constraint parameters. The † † \dagger symbol indicates that the corresponding value in the table is determined by equation ( 16 ).

				Inflexible Gen.			Intermittent Gen.						Flexible Gen.				Load Buses
Parameters		Units		1	7		2	3	4	5	6		8	9	10		11	12	13	14
${\alpha }_i$		$/MW		30	30		0	0	0	0	0		50	50	50		–	–	–	–
$\left\{d_i\right\}$		MW		0	21.7		87.8	0	9.0	3.5	13.5		94.2	11.2	0		7.6	29.5	6.1	14.9
$\mathrm{Im}\left\{d_i\right\}$		MVAR		0	12.7		-3.9	0	5.8	1.8	5.8		19	7.5	0		1.6	16.6	1.6	5
$\left\{g_i^{\max }\right\}$		MW		100	100		30	30	30	30	30		80	50	40		–	–	–	–
$\mathrm{Im}\left\{g_i^{\max }\right\}$		MVAR		10	50		$†$	$†$	$†$	$†$	$†$		40	24	24		–	–	–	–
$\left\{g_i^{\min }\right\}$		MW		0	0		0	0	0	0	0		0	0	0		–	–	–	–
$\mathrm{Im}\left\{g_i^{\min }\right\}$		MVAR		0	-40		$†$	$†$	$†$	$†$	$†$		0	-6	-6		–	–	–	–
$\left\{r_i^{\max }\right\}=-\left\{r_i^{\min }\right\}$		MW		0	0		$\infty$	$\infty$	$\infty$	$\infty$	$\infty$		$\infty$	$\infty$	$\infty$		–	–	–	–
$\mathrm{Im}\left\{r_i^{\max }\right\}=-\mathrm{Im}\left\{r_i^{\min }\right\}$		MVAR		0	0		$\infty$	$\infty$	$\infty$	$\infty$	$\infty$		$\infty$	$\infty$	$\infty$		–	–	–	–
Bus index (per [34])				1	2		4	7	10	11	13		3	6	8		5	9	12	14

Var/Symb		Space		Description
$n$		$𝐍$		number of network buses
$m$		$𝐍$		number of constraints of RAC-OPF
$k$		$𝐍$		ambient dimension of uncertainty set $\Xi$
$\ell$		$𝐍$		number of ellipsoids describing $\Xi$
$𝒢=(𝒱,ℰ)$		–		graph describing the power network
$𝒱$		–		set of vertices indexing the transmission buses
$ℰ$		–		set of edges indexing the transmission lines
$Y$		${𝐂}^{n\times n}$		bus admittance matrix
$𝝃$		–		random vector
$\xi$		${𝐑}^k$		a realization of the random vector $𝝃$
$\Xi$		–		support of the random vector $𝝃$
$\mu$		${𝐑}^k$		first-order moment of $𝝃$
$M$		${𝐑}^{k\times k}$		second-order moment of $𝝃$
${\ell }_{ij}^{\max }$		$𝐑$		active power flow capacity of transmission line $(i,j)$
${\alpha }_i$		${𝐑}_{+}$		marginal cost of active power generation at bus $i$
$\beta$		${𝐑}_{+}$		value of lost load (VOLL)
$v$		${ℒ}_{k,n}^2$		real-time bus voltage phasor
$g_i^0$		$𝐂$		day-ahead dispatch of generator $i$
$g_i$		${ℒ}_{k,1}^2$		real-time dispatch of generator $i$
$g_i^{\max }$		$𝐂$		nominal maximum capacity of generator $i$
$g_i^{\min }$		$𝐂$		nominal minimum capacity of generator $i$
${\overline{g}}_i$		${ℒ}_{k,1}^2$		real-time maximum capacity of generator $i$
${\underset{¯}{g}}_i$		${ℒ}_{k,1}^2$		real-time minimum capacity of generator $i$
$d_i$		$𝐂$		real-time demand at bus $i$
$r_i^{\min }$		$𝐂$		ramp-down limit of generator $i$
$r_i^{\max }$		$𝐂$		ramp-up limit of generator $i$
${\gamma }_i^{\max }$		$𝐂$		guaranteed maximum available power supply at bus $i$
${\gamma }_i^{\min }$		$𝐂$		guaranteed minimum available power supply at bus $i$

Equations

$$s_i=v^H{S}_{i}v,$$

$$s_{ij}=v^H{S}_{ij}v,$$

$$v_i^{\min }\le |v_i|\le v_i^{\max },$$

$$-{\ell }_{ij}^{\max }\le v^H{P}_{ij}v\le {\ell }_{ij}^{\max }.$$

$$\mu :=\mathbb{E}[\mathbf{\xi }]\text{and}M:=\mathbb{E}[\mathbf{\xi }{\mathbf{\xi }}^T],$$

$$\Xi :=\{\xi \in {\mathbf{R}}^k|{\xi }_1=1,{\xi }^T{W}_j\xi \ge 0,j\in [\ell ]\}.$$

$$W_j:=\left[\begin{array}{cc}{\omega }_j& w_j^T\\ w_j& -{\Omega }_j^T{\Omega }_j\end{array}\right],$$

$${\alpha }_i\Re \{g_i(\mathbf{\xi })\},i\in \mathcal{V},$$

$${\underline{g}}_i(\mathbf{\xi })\le g_i(\mathbf{\xi })\le {\overline{g}}_i(\mathbf{\xi }),i\in \mathcal{V}.$$

$$g_i^{\min }\le {\underline{g}}_i(\mathbf{\xi })\le {\overline{g}}_i(\mathbf{\xi })\le g_i^{\max },i\in \mathcal{V}.$$

$$r_i^{\min }\le g_i(\mathbf{\xi })-g_i^0\le r_i^{\max },i\in \mathcal{V},$$

$$g_i^{\min }\le g_i^0\le g_i^{\max },i\in \mathcal{V}.$$

$$g_i(\mathbf{\xi })=g_i^0.$$

$$g_i^{\min }\le g_i(\mathbf{\xi })\le g_i^{\max }.$$

$${\underline{g}}_i(\mathbf{\xi })\le g_i(\mathbf{\xi })\le {\overline{g}}_i(\mathbf{\xi }).$$

$$\mathrm{minimize}$$

$$g^0\in {\mathbf{C}}^n,g\in {\mathcal{L}}_{k,n}^2,v\in {\mathcal{L}}_{k,n}^2$$

$$g_i^{\min }\le g_i^0\le g_i^{\max },$$

$${\underline{g}}_i(\xi )\le g_i(\xi )\le {\overline{g}}_i(\xi ),$$

$$r_i^{\min }\le g_i(\xi )-g_i^0\le r_i^{\max },$$

$$g_i(\xi )-v(\xi {)}^H{S}_{i}v(\xi )=d_i,$$

$$v_i^{\min }\le |v_i(\xi )|\le v_i^{\max },$$

$$|v(\xi {)}^H{P}_{ij}v(\xi )|\le {\ell }_{ij}^{\max },$$

$$\displaystyle\max_{\xi\in\Xi}\ \Bigg{\{}\sum_{i=1}^{n}\alpha_{i}{\rm Re}\{g_{i}(\xi)\}\Bigg{\}}.$$

$$\displaystyle\underset{}{\text{minimize}}\ \ \$$

$$x\in {\mathbf{R}}^{2n},v\in {\mathcal{L}}_{k,n}^2$$

$$v(\xi {)}^H{A}_{i}v(\xi )+b_i^{T}x\le c_i^T\xi ,i\in [m],\forall \xi \in \Xi$$

$$Ex\le f,$$

$$v(\mathbf{\xi })=V\mathbf{\xi },$$

$$\displaystyle\underset{}{\text{minimize}}\ \ \$$

$$x\in {\mathbf{R}}^{2n},V\in {\mathbf{C}}^{n\times k}$$

$${\xi }^T{V}^H{A}_{i}V\xi +b_i^{T}x\le c_i^T\xi ,i\in [m],\forall \xi \in \Xi$$

$$Ex\le f,$$

$$\mathbb{E}[{\mathbf{\xi }}^T{V}^H{A}_{0}V\mathbf{\xi }]=\mathbb{E}[\mathrm{tr}(\mathbf{\xi }{\mathbf{\xi }}^T{V}^H{A}_{0}V)]=\mathrm{tr}(\mathbb{E}[\mathbf{\xi }{\mathbf{\xi }}^T]V^H{A}_{0}V).$$

$$\mathrm{minimize}$$

$$x\in {\mathbf{R}}^{2n},V\in {\mathbf{C}}^{n\times k},\Lambda \in {\mathbf{R}}^{m\times \ell }$$

$$V^H{A}_{i}V-C_i+(b_i^{T}x)e_1{e}_1^T-\sum _{j=1}^{\ell }[\Lambda {]}_{ij}{W}_j$$

$$\forall i\in [m],$$

Full text

I Introduction

The AC optimal power flow (AC-OPF) problem is a fundamental decision problem that is at the heart of power system operations [2]. The AC-OPF problem is a nonconvex optimization problem, where the objective is to minimize the cost of generation subject to power balance constraints described by Kirchhoff’s current and voltage laws, and operational constraints reflecting real and reactive limits on power generation, bus voltage magnitudes, and power flows along transmission lines flows. It is also common to enforce contingency constraints to ensure that the power system can withstand sudden disturbances, such as generator or line outages [3]. The nonconvexity of the AC-OPF problem is in part due to the need to enforce quadratic constraints, which are indefinite in the vector of complex bus voltages. The treatment of such nonconvexities in the AC-OPF problem has traditionally relied on the use of local constrained-optimization methods, or the use of approximate linear models of power flow (e.g., DC-OPF) to convexify the feasible region of the underlying optimization problem[4]. More recently, considerable effort has been made to identify conditions under which an optimal solution to AC-OPF can be obtained from a solution to its semidefinite programming relaxation [5, 6, 7, 4]. According to the US Federal Energy Regulatory Commission, a 5% increase in the efficiency of algorithms for AC-OPF will yield six billion dollars in savings per year in the United States alone [8].

Recently, growing concerns about climate change have led many US states to implement policies, which require that a large fraction of their electricity come from renewable energy resources such as wind and solar. A fundamental challenge facing the deep integration of such resources stems from the need to accommodate the intrinsic uncertainty in their power supply. Doing so efficiently will require the development of robust optimization methods for the AC optimal power flow (AC-OPF) problem. In its most elemental formulation, the robust AC optimal power flow (RAC-OPF) problem gives rise to a two-stage robust optimization problem, in which the system operator must determine a day-ahead generation schedule that minimizes the expected cost of dispatch, given an opportunity for recourse to adjust its day-ahead schedule in real-time when the uncertain system variables have been realized (e.g., the power that can be supplied from wind and solar resources). The RAC-OPF problem is a nonconvex, infinite-dimensional optimization problem in its most general form due to the nonconvexity of the underlying AC power flow constraints, and the need to optimize over an infinite-dimensional recourse policy space.

Related Work:  In order to treat the nonconvexity of the RAC-OPF problem, a large fraction of the literature prescribes techniques that rely on a DC linear approximation of the power flow model [9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. The computational intractability associated with the need to optimize over an infinite-dimensional recourse policy space is predominantly addressed by employing affine or piecewise-affine approximations of the recourse policy space [19, 9, 10, 11, 12, 13, 14, 15, 16, 17, 20, 18]. Lastly, both robust [19, 10, 11, 14, 21, 18] and chance-constrained [9, 17, 22, 15, 16, 20] formulations have been proposed to treat the uncertainty in the constraints, which define the RAC-OPF problem. As the title of this paper suggests, we adopt a robust approach to constraint satisfaction.

The chance-constrained paradigm assigns a distribution to the uncertain variables, and the system constraints are enforced up to a prespecified probability level. The procurement of such distributions, however, is challenging as distributions describing uncertain power system parameters, such as renewable energy generation, can be difficult to identify [23]. In order to account for the potential inaccuracy in the specification of the underlying distribution, the authors of [24] consider ambiguous chance constraint formulations using DC power flow models. Under this paradigm, the underlying distribution is assumed to belong to a closed ball centered around a known distribution. An alternative treatment of ambiguous chance constraints is provided in [25]. In addition to the basic difficulty inherent to the identification of an accurate distributional model, chance constraints do not explicitly account for the magnitude of constraint violations when they occur—however small their a priori specified probability. To account for this, the conditional value at risk (CVaR) is commonly used in the place of chance constraints to promote solutions that minimize the expected magnitude of such constraint violations [17, 16]. In contrast to the chance-constrained paradigm, robust optimization takes a deterministic approach to uncertainty modeling. In particular, the uncertain parameters are assumed to vary in a known and bounded uncertainty set. A robust solution is one that optimizes the objective function and remains feasible for any realization of the uncertain parameters in the given set. This immunity of robust solutions, however, comes at the expense of potential conservatism in their performance. Under certain mild assumptions, robust linear constraints arising in approximate DC power flow models admit equivalent reformulations as finite-dimensional conic linear constraints [10]. In general, however, robust nonlinear constraints, which arise in AC power flow models, do not admit such equivalent reformulations. The predominant approach to their treatment relies on scenario or sample-based approximation techniques [26], which give rise to outer approximations (relaxations) of the robustly feasible set.

Summary of Results:   In this paper, we formulate RAC-OPF as a two-stage robust optimization problem with recourse. The formulation considered departs from the majority of the extant literature given its treatment of the full AC power flow model. To our knowledge, the only other papers in the literature that treat the AC model in the robust optimization framework are [19, 27, 20, 21]. They address the nonconvexity, which arises from the AC power flow equations, by means of a convex (second-order cone or semidefinite) relaxation.111We refer the reader to [6, 7, 4] for a detailed derivation and analysis of semidefinite and second-order cone relaxations for the deterministic (zero-recourse) AC-OPF problem. A crucial assumption made in these papers is the exactness of their convex relaxations. Exactness of such relaxations for the RAC-OPF problem is not guaranteed, and, in particular, the solutions generated by these relaxations are not guaranteed to be feasible for RAC-OPF. In this paper, we adopt an approach that relies on the restriction of the space of recourse policies to those which are affine in the uncertain problem data. Under this restriction and the assumption that there is adjustable generation or load at every bus in the power transmission network, we develop a technique to approximate the RAC-OPF problem from within by a finite-dimensional semidefinite program. We establish a sufficient condition under which the resulting semidefinite program—a convex inner approximation to RAC-OPF—yields recourse policies, which are guaranteed to be feasible for RAC-OPF.

Organization:  The paper is organized as follows. In Sections II and III, we develop the power system model and provide a detailed formulation of RAC-OPF, respectively. In Section IV, we offer a detailed derivation of the semidefinite programming inner approximation of RAC-OPF, and provide a sufficient condition under which the resulting approximation is guaranteed to have a nonempty feasible region. In Section V, we describe an iterative optimization method that generates a sequence of feasible affine recourse policies with nonincreasing costs. Finally, we illustrate the effectiveness of the proposed optimization method on the WSCC 9-bus and IEEE 14-bus power systems with different levels of renewable resource penetration and uncertainty. Section VII concludes the paper. Detailed proofs for the theoretical results contained in this paper can be found in the Appendix.

Notation:  Let $\mathbf{R}$, $\mathbf{C}$, and $\mathbf{N}$ denote the set of real, complex, and natural numbers, respectively. Given $m\in\mathbf{N}$, define the set of natural numbers from 1 to $m$ by $[m]\coloneqq\{1,2,\dots,m\}$. Let $e_{i}$ denote the real $i^{\rm th}$ standard basis vector, whose dimension will be apparent from the context. Given any pair of complex numbers $z_{1},z_{2}\in\mathbf{C}$, we write $z_{1}\leq z_{2}$ if and only if $\real\{z_{1}\}\leq\real\{z_{2}\}\text{ and }{\rm Im}\{z_{1}\}\leq{\rm Im}\{z_{2}\}$. Given a matrix $X\in\mathbf{C}^{m\times n}$, we denote its conjugate transpose by $X^{*}$, and its $(i,j)$ entry by $[X]_{ij}$. Given a matrix $X\in\mathbf{R}^{m\times n}$, we denote its transpose by $X^{T}$ and use $X\geq 0$ to mean that $X$ is entrywise nonnegative. We denote the trace of a matrix $X$ by ${\rm tr}(X)$. Denote by $I_{n}$ the $n\times n$ identity matrix. Denote by $\mathbf{H}^{n}$ the set of $n\times n$ Hermitian matrices. We use $X\succeq 0$ to mean that the matrix $X\in\mathbf{H}^{n}$ is positive semidefinite. Finally, let $\mathcal{L}_{k,n}^{2}$ denote the space of all square-integrable, Borel measurable functions from $\mathbf{R}^{k}$ to $\mathbf{C}^{n}.$ A complex-valued function $f$ on $\mathbf{R}^{k}$ is said to be Borel measurable if both $\real\{f\}$ and ${\rm Im}\{f\}$ are real-valued Borel measurable. We summarize frequently used symbols and variables in Appendix E.

II Power System Model

In this Section, we develop the robust AC optimal power flow (RAC-OPF) problem. We consider a power system consisting of a heterogeneous mix of generators, which differ in terms of their predictability and controllability. We adopt the perspective of the independent system operator (ISO), whose aim is to dispatch generators in order to minimize the expected cost of serving demand, while respecting generation and transmission capacity constraints—an optimization problem that belongs to the class of security-constrained optimal power flow problems. We consider an optimization model that consists of two stages: day-ahead (DA) and real-time (RT). In the day-ahead stage, the ISO determines an initial dispatch of generators subject to uncertainty in certain system variables, e.g., the supply that will be available from renewable energy resources in real-time. Day-ahead scheduling decisions are critical, as certain generators (e.g., nuclear and coal) have ramping constraints that limit the extent to which they can adjust their power output in real-time. Accordingly, such ramp-constrained generators must be scheduled to produce well in advance of the delivery time, and, therefore, prior to the realization of certain a priori uncertain variables, e.g., wind power availability. In the real-time stage, all a priori uncertain variables are realized, and the ISO is given a recourse opportunity to adjust its DA schedule to balance the system at minimum cost. The ramping constraints dictate the extent to which each generator can adjust its power output around its day-ahead set-point in real-time. Essentially, the calculation of a DA schedule that minimizes the expected cost of generation—subject to optimal recourse in real-time—entails the solution of a robust optimization problem with recourse. We provide a precise formulation of this problem in Section III.

II-A AC Power Flow Model

We consider a power transmission network described by an undirected graph ${\cal G}\coloneqq({\cal V},{\cal E})$. The set of vertices ${\cal V}:=[n]$ index the transmission buses, and the set of edges ${\cal E}\subseteq{\cal V}\times{\cal V}$ index the transmission lines between buses. We assume that $(j,i)\in{\cal E}$ if and only if $(i,j)\in{\cal E}$.

The AC power flow equations are formulated according to Kirchhoff’s voltage and current laws, which relate bus power injections to voltages [28]. We denote by $Y\in\mathbf{C}^{n\times n}$ the bus admittance matrix, $s\in\mathbf{C}^{n}$ the vector of complex (net) bus power injections (generation minus demand), and $v\in\mathbf{C}^{n}$ the vector of complex bus voltages. The AC power balance equations are given by

[TABLE]

where $S_{i}\coloneqq Y^{H}e_{i}e_{i}^{T}$ for all $i\in{\cal V}$. For each line $(i,j)\in{\cal E}$, we denote by $s_{ij}\in\mathbf{C}$ the complex power flow from bus $i$ to bus $j$. It satisfies

[TABLE]

where $S_{ij}\coloneqq e_{i}e_{i}^{T}(\widehat{y}_{ij}/2-[Y]_{ij})^{H}+e_{j}e_{i}^{T}[Y]_{ij}^{H}$ for all $(i,j)\in{\cal E}$. Here, $\widehat{y}_{ij}\in\mathbf{C}$ represents the total shunt admittance of line $(i,j)$.

We require that the following constraints be enforced. The first set of constraints limit the range of acceptable voltage magnitudes at each bus $i\in{\cal V}$,

[TABLE]

where $v_{i}^{\min}\in\mathbf{R}$ and $v_{i}^{\max}\in\mathbf{R}$ denote lower and upper limits, respectively, on the bus $i$ voltage magnitude. We also consider active power flow capacity constraints on the transmission lines. Namely, for each transmission line $(i,j)\in{\cal E}$, the active power flow from bus $i$ to bus $j$ is required to satisfy

[TABLE]

Here, we define $P_{ij}\coloneqq(S_{ij}+S_{ij}^{H})/2$, and let $\ell_{ij}^{\max}\in\mathbf{R}$ denote the active power flow capacity of transmission line $(i,j)$. We note that—within the framework of this paper—it is also possible to accommodate transmission line constraints that limit the magnitude of the current flowing through a line, as these also give rise to quadratic inequality constraints in the vector of complex bus voltages. The treatment of apparent power flow constraints, however, is beyond the scope of the present paper, as these amount to quartic inequality constraints in the vector of complex bus voltages.

II-B Uncertainty Model

All of the ‘uncertain’ quantities appearing in this paper are described according to the random vector $\bm{\xi}$, which is defined according to the probability space $(\mathbf{R}^{k},{\cal B}(\mathbf{R}^{k}),\mathbb{P})$. Here, the Borel $\sigma$-algebra ${\cal B}(\mathbf{R}^{k})$ is the set of all events that are assigned probabilities by the measure $\mathbb{P}$. We denote the first and second-order moments of $\bm{\xi}$ by

[TABLE]

where $\mathbb{E}[\cdot]$ denotes the expectation operator with respect to $\mathbb{P}$. Adopting a standard notational convention, we will use $\xi$ (normal face) to denote realizations taken by the random vector $\bm{\xi}$ (bold face). We assume throughout the paper that the support of the random vector $\bm{\xi}$ is nonempty, compact, and representable as

[TABLE]

Here, each matrix $W_{j}\in\mathbf{R}^{k\times k}$ is defined according to

[TABLE]

where $\omega_{j}\in\mathbf{R}$, $w_{j}\in\mathbf{R}^{k-1}$, and $\Omega_{j}\in\mathbf{R}^{n_{j}\times(k-1)}$ for some $n_{j}\in\mathbf{N}$. It is important to note that the representation of the support in (3) is general enough to describe any subset of the hyperplane $\{\xi\in\mathbf{R}^{k}\ |\ \xi_{1}=1\}$ that is defined according to a finite intersection of arbitrary half spaces and ellipsoids. We will occasionally refer to the support set $\Xi$ as the ‘uncertainty set’ associated with the random vector $\bm{\xi}$. Naturally, the particular specification of the uncertainty set will depend on the specific nature of the uncertain quantities that the random vector $\bm{\xi}$ is meant to capture, e.g., the maximum active power that a wind generator is able to produce at a node.

Some remarks regarding our uncertainty model are in order. First, the requirement that $\xi_{1}=1$ for all $\xi\in\Xi$ is for notational convenience, as it allows one to represent affine functions of $(\xi_{2},\dots,\xi_{k})$ as linear functions of $\xi$. Second, it is important to emphasize that all of the results contained in this paper depend on the probability distribution of the random vector $\bm{\xi}$ only through its support, mean, and second-order moment. No additional information about the distribution is required. The following mild technical assumption is assumed to hold throughout the paper.

Assumption 1**.**

There exists $\xi\in\Xi$ such that $\xi^{T}W_{j}\xi>0$ for all $j\in[\ell]$.

The assumption that the support set $\Xi$ admits a strictly feasible point will prove useful to the derivation of our subsequent theoretical results, as it ensures that $\Xi$ spans all of $\mathbf{R}^{k}$. This, in turn, guarantees that the second-order moment matrix $M$ is positive definite and invertible. We refer the reader to [29, Prop. 2] for a proof of this claim.

II-C Generator and Load Models

Load Model:   The real-time demand for power at each bus $i\in{\cal V}$ is assumed to be fixed and known. We denote it by $d_{i}\in\mathbf{C}$.

Generator Model:  For ease of exposition, we assume throughout the paper that there is at most a single generator at each bus $i\in{\cal V}$. We consider a generator model in which the real-time (RT) supply of power, as determined by the ISO, is allowed to depend on the realization of the random vector $\bm{\xi}$. Hence, we denote the power produced by generator $i$ in real-time by $g_{i}(\bm{\xi})$, where $g_{i}\in\mathcal{L}^{2}_{k,1}$ is a recourse function determined by the ISO for each generator $i$. The cost incurred by each generator $i$ for producing $g_{i}(\bm{\xi})$ is assumed to be linear in the active power produced. It is defined as

[TABLE]

where $\alpha_{i}\geq 0$ denotes the marginal cost of real power generation at bus $i\in{\cal V}$.

We consider a generator model in which the power capacity available to each generator in real-time is allowed to be uncertain day-ahead. Namely, the power produced by each generator $i\in{\cal V}$ in real-time must satisfy

[TABLE]

Here, $\underline{g}_{i}\in\mathcal{L}^{2}_{k,1}$ and $\overline{g}_{i}\in\mathcal{L}^{2}_{k,1}$ denote lower and upper bounds on the power produced by generator $i$ in real-time. Such a model is general enough to describe apriori uncertainty in renewable power supply, as well as unscheduled generator outages. We assume that the random generation capacities satisfy

[TABLE]

Here, $g_{i}^{\min}\in\mathbf{C}$ and $\ g_{i}^{\max}\in\mathbf{C}$ denote the nameplate minimum and maximum capacities of generator $i$, respectively. The corresponding vectors are denoted by $\underline{g}(\bm{\xi})$, $\overline{g}(\bm{\xi})$, $g^{\min}$, and $g^{\max}$.

In practice, generators have ramping constraints that limit the extent to which they can adjust their power production in real-time. We model the limited ramping capability of each generator $i$ according to the following pair of constraints

[TABLE]

where $r^{\min}_{i}\in\mathbf{C}$ and $r^{\max}_{i}\in\mathbf{C}$ represent the ramp-down and ramp-up limits, respectively, associated with generator $i$. Here, $g_{i}^{0}\in\mathbf{C}$ denotes generator $i$’s day-ahead (DA) dispatch, also determined by the ISO. The DA dispatch of each generator is required to satisfy its nameplate generation capacity constraints given by

[TABLE]

Example 1** (Generator types).**

The generator model that we consider in this paper captures a wide range of generator types. We provide several important examples below. Let $g_{i}^{0}$ be a DA dispatch level satisfying (7). Generator $i$ is said to be:

• •

Completely inflexible (e.g., nuclear) if its RT power output is restricted to

[TABLE]

• •

Completely flexible (e.g., gas, oil) if its RT power output is restricted to

[TABLE]

• •

Intermittent (e.g., wind, solar) if its RT power output is restricted to

[TABLE]

We state the following technical assumption, which requires that the RT generation capacities exhibit a linear dependence on the random vector $\bm{\xi}$. Assumption 2 is required to hold for the remainder of the paper.

Assumption 2**.**

There exist matrices $\underline{G}\in\mathbf{C}^{n\times k}$ and $\overline{G}\in\mathbf{C}^{n\times k}$ such that $\underline{g}(\bm{\xi})=\underline{G}\bm{\xi}$ and $\overline{g}(\bm{\xi})=\overline{G}\bm{\xi}$.

III Robust AC Optimal Power Flow

Building on the previously defined models, we formulate the robust AC optimal power flow (RAC-OPF) problem as follows.

[TABLE]

As previously described, the RAC-OPF problem amounts to a two-stage robust optimization problem with recourse. The single-period formulation of RAC-OPF that we consider is similar in structure to the single-period formulations studied in [19, 10, 30, 20]. We briefly summarize the timing and structure of the decision variables and constraints of the RAC-OPF problem.

• •

The first-stage (day-ahead) decisions entail the determination of a DA generator dispatch $g_{0}\in\mathbf{C}^{n}$ subject to optimal recourse in the second stage, which will adjust the DA dispatch given a realization of the random vector $\bm{\xi}$.

• •

In the second-stage (real-time), the random vector $\bm{\xi}$ is realized, and the ISO is given a recourse opportunity to adjust its DA generator dispatch to balance the system at minimum cost. The second-stage decision entails the determination of the RT generator dispatch $g\in\mathcal{L}^{2}_{k,n}$ and the RT bus voltages $v\in\mathcal{L}^{2}_{k,n}$.

• •

All decisions must be jointly determined in such a manner as to (i) minimize the expected cost of generation, and (ii) guarantee that all system constraints are satisfied given any realization $\xi\in\Xi$ of the random vector $\bm{\xi}$ in real-time, i.e., robust constraint satisfaction.

Remark 1** (Minimax formulation).**

While the formulation of the RAC-OPF problem in (8) entails minimizing the expected cost of generation, the computational methods and theoretical results developed in Sections IV -V can be generalized to accommodate a minimax formulation of the RAC-OPF problem, which entails minimizing the maximum (worst-case) cost of generation:

[TABLE]

The resulting minimax formulation of the RAC-OPF problem under the objective function (9) can be equivalently reformulated to resemble (8) by putting it in its epigraph form.

III-A Concise Formulation of RAC-OPF

It will be convenient to our analysis in the sequel to work with a more concise representation of the RAC-OPF problem. We do so by first eliminating the RT generator dispatch variables $g\in\mathcal{L}^{2}_{k,n}$ through their direct substitution according to the nodal power balance equations. Second, by redefining the DA generator dispatch $g_{0}\in\mathbf{C}^{n}$ as a real vector $x:=[\real\{g_{0}\}^{T},{\rm Im}\{g_{0}\}^{T}]^{T}$, one can rewrite problem (8) more compactly in the following form:

[TABLE]

where $m:=10n+2|{\cal E}|$. We remark that in the above reformulation of problem (8), we have eliminated the constant term $\sum_{i=1}^{n}\alpha_{i}\real\{d_{i}\}$ from the objective function as this does not affect the optimal solution of the RAC-OPF problem. It is straightforward to construct the matrices $E\in\mathbf{R}^{4n\times 2n}$, $f\in\mathbf{R}^{4n}$, $A_{i}\in\mathbf{H}^{n}$ ($i=0,\dots,m$), $b_{i}\in\mathbf{R}^{2n}$ ($i=1,\dots,m$), and $c_{i}\in\mathbf{R}^{k}$ ($i=1,\dots,m$) given the underlying problem data specified in the RAC-OPF problem (8). We refer the reader to Appendix A for their specification.

Remark 2** (Eliminating quadratic equality constraints).**

We remark that the formulation of the original RAC-OPF problem (8) assumes that there is adjustable generation at every bus in the power transmission network. The advantage of this rather limiting assumption is that it enables the elimination of all nodal power balance equality constraints in the equivalent reformulation of the RAC-OPF problem given by ${\cal P}$. The ability to eliminate these nonconvex quadratic equality constraints will be essential to the convex inner approximation technique developed in Section IV-B. In Appendix D, we provide an alternative formulation of the RAC-OPF problem to accommodate the treatment of more general power systems in which load shedding is permitted at non-generator (load) buses, where any reduction in load is penalized according to the value of lost load (VOLL)). We refer the reader to Section VI for several numerical case studies, which asses the extent to which the allowance of load shedding manifests in an actual reduction in load under the dispatch policies proposed in this paper. Finally, we note that the approximation technique proposed in this paper cannot be applied to power systems with transmission buses that have neither adjustable generation or load.

IV Convex Inner Approximation of RAC-OPF

Problem ${\cal P}$ is computationally intractable, in general, as it is both infinite-dimensional and nonconvex. The nonconvexity is due, in part, to the feasible set, which is defined by a number of indefinite quadratic inequality constraints in the vector of complex bus voltages. The infinite-dimensionality of the optimization problem ${\cal P}$ derives from both the infinite-dimensionality of the recourse decision variables, and the infinite number of constraints due to the infinite cardinality of the uncertainty set $\Xi$. In what follows, we develop a systematic approach to approximate problem ${\cal P}$ from within by a finite-dimensional semidefinite program, and provide a sufficient condition under which the resulting inner approximation is guaranteed to have a nonempty feasible region. The proposed method for approximation centers on the restriction of the infinite-dimensional space of recourse policies to those which are linear in the random vector $\bm{\xi}$.

IV-A Affine Recourse Policies

As the initial step in the derivation of a tractable inner approximation to problem ${\cal P}$, we first restrict the functional form of the recourse decision variables (i.e., the complex bus voltages) to be linear in the random vector $\bm{\xi}$.222We note that this restriction on the functional form of the complex bus voltages implies a quadratic dependency of the nodal power generation levels on the random vector $\bm{\xi}$. That is to say, $g_{i}(\bm{\xi})=\bm{\xi}^{T}V^{H}S_{i}V\bm{\xi}+d_{i}$ for each node $i\in{\cal V}$. That is to say, we require that

[TABLE]

where $V\in\mathbf{C}^{n\times k}$. This restriction to affine recourse policies gives rise to the following optimization problem IV-A, which stands as an inner approximation to the original problem ${\cal P}$.

[TABLE]

We have used linearity of the expectation and trace operators, and the invariance of trace under cyclic permutations to massage the original objective function to obtain

[TABLE]

The resulting problem IV-A amounts to a semi-infinite, nonconvex quadratically constrained quadratic program.333A semi-infinite program is an optimization problem with an infinite number of constraints, and finitely many decision variables. More specifically, the restriction to affine recourse policies yields an optimization problem that has finite-dimensional decision variables. However, Problem IV-A remains to be computationally intractable, as it requires the satisfaction of infinitely many constraints due to the continuous structure of the uncertainty set $\Xi$. We address this issue in Lemma 1 by employing weak duality to obtain a sufficient set of finitely many constraints. Such an approximation of the infinite constraint set can also be derived through a direct application of the so-called $S$-procedure [31]. We state Lemma 1 without proof, as it follows directly from Proposition 6 in [29].

Lemma 1**.**

Let $P\in\mathbf{H}^{k},\ q\in\mathbf{R}^{k}$, $r\in\mathbf{R}$, and $Q\coloneqq(e_{1}q^{T}+qe_{1}^{T})/2$. Consider the following two statements:

1. (i)

$\xi^{T}P\xi+q^{T}\xi+r\leq 0$* for all $\xi\in\Xi$,* 2. (ii)

$\exists$* $\lambda\in\mathbf{R}^{\ell}$ with $\lambda\leq 0$ and $P+Q+re_{1}e_{1}^{T}-\sum\limits_{j=1}^{\ell}\lambda_{j}W_{j}\preceq 0,$*

where $W_{j}$ is as defined in (4). For any $\ell\in\mathbf{N}$, it holds that (ii) implies (i). If $\ell=1$, then (i) and (ii) are equivalent.

Using Lemma 1, one can approximate the infinite constraint set of problem IV-A from within by finitely many matrix inequality constraints. More precisely, a direct application of Lemma 1 to each of the quadratic constraints in problem IV-A gives rise to the following finite-dimensional optimization problem:

[TABLE]

where we define $C_{i}\coloneqq(e_{1}c_{i}^{T}+c_{i}e_{1}^{T})/2$ for each $i\in[m]$.

Remark 3**.**

We remark that $\Lambda\in\mathbf{R}^{m\times\ell}$ is a decision variable in problem ${\cal P}_{\rm II}$, which emerges from the from the application of Lemma 1 to the robust inequality constraints in problem ${\cal P}_{I}.$ It follows from Lemma 1 that problem ${\cal P}_{\rm II}$ is an inner approximation to problem IV-A in general; and is equivalent to problem IV-A when $\ell=1$.

IV-B Convexifying the Inner Approximation

Problem ${\cal P}_{\rm II}$ is a finite-dimensional inner approximation to the original problem ${\cal P}$. It, however, remains to be nonconvex, because of the indefinite quadratic functions appearing in both the objective and the inequality constraints. In what follows, we develop a method to convexify problem ${\cal P}_{\rm II}$ from within by replacing each indefinite quadratic function with a majorizing convex quadratic function. We state the resulting convex program, which approximates ${\cal P}_{\rm II}$ from within, in Proposition 1.

The proposed method is based on the decomposition of an indefinite quadratic function as the difference of convex functions, i.e., the sum of a convex quadratic function and a concave quadratic function. We construct a convex global overestimator of the original indefinite quadratic function by linearizing the concave function at a point.444In Section IV-C, we provide a method to select the point around which the linearization is constructed in order to guarantee that the resulting convex feasible set is nonempty. More precisely, for each matrix $A_{i}$, define the decomposition

[TABLE]

where $A^{+}_{i}\succeq 0$ and $A^{-}_{i}\preceq 0$ denote the positive semidefinite and negative semidefinite parts of $A_{i}$, respectively. Using this matrix decomposition, define the function $H_{i}:\mathbf{C}^{n\times k}\times\mathbf{C}^{n\times k}\rightarrow\mathbf{H}^{k}$ according to

[TABLE]

for each $i\in[m]$. The first term of $H_{i}$ is the convex component of the original quadratic function $V^{H}A_{i}V$. The remaining terms represent the linearization of the concave component at a point $Z$. Consequently, for any matrix $Z$, the function $H_{i}(V,Z)$ is matrix convex in $V$.555A function $f:\mathbf{C}^{n\times k}\rightarrow\mathbf{H}^{k}$ is said to be matrix convex if for all matrices $X,Y$ and $0\leq\theta\leq 1$, we have $f(\theta X+(1-\theta)Y)\preceq\theta f(X)+(1-\theta)f(Y)$. The following result highlights two important properties of $H_{i}$. Its proof can be found in Appendix B-A.

Lemma 2**.**

Let $Z\in\mathbf{C}^{n\times k}$. For each $i\in[m]$, it holds that

1. (i)

$V^{H}A_{i}V\preceq H_{i}(V,Z),\quad\forall\ V\in\mathbf{C}^{n\times k}$, 2. (ii)

${\rm tr}(MV^{H}A_{i}V)\leq{\rm tr}(MH_{i}(V,Z)),\quad\forall\ V\in\mathbf{C}^{n\times k}$.

Property (i) provides a way of approximating the nonconvex feasible set of problem ${\cal P}_{\rm II}$ from withing by a convex set. Property (ii), on the other hand, provides way of majorizing the nonconvex objective of problem ${\cal P}_{\rm II}$ with a convex function. In Proposition 1, we employ these approximations to specify a convex program whose optimal solution is guaranteed to be a feasible solution for the original problem ${\cal P}$. Its proof follows directly from Lemma 2. We, therefore, omit it for the sake of brevity.

Proposition 1**.**

Let $V_{0}\in\mathbf{C}^{n\times k}$, and suppose that $(\overline{x},\overline{V},\overline{\Lambda})$ is an optimal solution for the following convex program:

[TABLE]

Define the function $\overline{v}\in{\cal L}^{2}_{k,n}$ according to $\overline{v}(\bm{\xi})=\overline{V}\bm{\xi}$. Then $(\overline{x},\overline{v})$ is a feasible solution for the original problem ${\cal P}$.

We note that problem ${\cal P}_{\rm III}(V_{0})$ can be equivalently reformulated as a semidefinite program using Schur’s complement formula. We refer the reader to Appendix C for the details of this reformulation.

IV-C Guaranteeing Nonemptiness of the Inner Approximation

In order to convexify the RAC-OPF problem according the method developed in Section IV-B, one has to select a matrix $V_{0}\in\mathbf{C}^{n\times k}$, which yields an inner approximation ${\cal P}_{\rm III}(V_{0})$ with a nonempty feasible set. In what follows, we develop a method to compute such a matrix.

The method we propose entails the calculation of a day-ahead dispatch $g^{0}\in\mathbf{C}^{n}$, which is guaranteed to be feasible for the RAC-OPF problem without requiring adjustment (recourse) in real-time. In order to do so, we must first characterize the guaranteed range of available power supply at each bus. For each bus $i\in{\cal V}$, this amounts to the specification of upper and lower limits $\gamma_{i}^{\max}\in\mathbf{C}$ and $\gamma_{i}^{\min}\in\mathbf{C}$, such that

[TABLE]

We specify these limits according to

[TABLE]

It is important to note that the limits $\gamma_{i}^{\min}$ and $\gamma_{i}^{\max}$ can be efficiently calculated, as the optimization problems in (11) are convex quadratically constrained quadratic programs. This follows from the assumed linearity of the objective function in $\xi$ (cf. Assumption 2), and the definition of $\Xi$ as the finite intersection of ellipsoids and half spaces (cf. Eq. (3)).

Using these conservative generation limits, a day-ahead dispatch that is guaranteed to be feasible for the RAC-OPF problem can calculated by solving the following (deterministic) zero-recourse AC-OPF problem.

[TABLE]

The following result shows that any feasible solution to the zero-recourse AC-OPF problem (12) can be used to construct a matrix $V_{0}$ that is guaranteed to induce a convex inner approximation ${\cal P}_{\rm III}(V_{0})$ of the RAC-OPF problem with a nonempty feasible region. The proof of Proposition 2 can be found in Appendix B-B.

Proposition 2**.**

Let $v_{0}\in\mathbf{C}^{n}$ be a feasible solution to (12), and define a matrix $V_{0}\coloneqq v_{0}e_{1}^{T}$. It follows that the optimization problem ${\cal P}_{\rm III}(V_{0})$ has a nonempty feasible region.

Several comments are in order. First, it is important to note that, despite being deterministic, the zero-recourse AC-OPF problem (12) is nonconvex and computationally intractable, in general. However, there are many off-the-shelf optimization routines (e.g., Matpower [33]) that are effective in producing feasible solutions to problem (12). Second, it is also crucial to emphasize that Proposition 2 is only useful if the zero-recourse AC-OPF problem (12) is in fact feasible. A necessary condition for the feasibility of problem (12) is that $\gamma_{i}^{\min}\leq\gamma_{i}^{\max}$ for every generator $i\in{\cal V}$. This condition is clearly satisfied by conventional generators with ‘firm’ real-time generation capacities, i.e., $\underline{g}_{i}(\xi)=g_{i}^{\min}$ and $\overline{g}_{i}(\xi)=g_{i}^{\max}$ for all $\xi\in\Xi$. This condition is also satisfied by intermittent generators whose real-time active power supply can be fully ‘curtailed’, i.e., $\real\{\underline{g}_{i}(\xi)\}=0$ for all $\xi\in\Xi$. Such an assumption of curtailable supply is reasonable for modern-day wind and solar power facilities. We refer the reader to Section VI-A for a specific example of curtailable intermittent generation. It is also important to note that the assumption that intermittent generators be curtailable is necessary, as problem (12) may become infeasible for power systems with intermittent generators whose real-time supply is uncertain in day-ahead and cannot be curtailed in real-time, i.e., $\underline{g}_{i}(\xi)=\overline{g}_{i}(\xi)$ for all $\xi\in\Xi$.

V Sequential Convex Approximation Method

In what follows, we describe a recursive method that builds upon our previous development to generate a sequence of cost-improving convex inner approximations to the RAC-OPF problem. Let $v_{0}\in\mathbf{C}^{n}$ be a feasible solution to the zero-recourse AC-OPF problem (12), and define the matrix $V_{0}\coloneqq v_{0}e_{1}^{T}$. Consider a recursion of the form

[TABLE]

Here, $\mathcal{F}(V_{t})$ denotes the feasible set of problem ${\cal P}_{\rm III}(V_{t})$, which is parameterized by the matrix $V_{t}$. The recursive algorithm (13) can be interpreted as a successive convex majorization-minimization method. An optimal solution $(x_{t},V_{t})$ associated with each step $t$ of the recursive method can be mapped to a feasible solution $(x_{t},v_{t}(\bm{\xi}))$ of the RAC-OPF problem (8), where $v_{t}(\bm{\xi})=V_{t}\bm{\xi}$ is an affine recourse policy in $\bm{\xi}$. Proposition 3 establishes two important properties of the recursive method. The recursive method is (i) guaranteed to yield a nonempty convex inner approximation to the RAC-OPF problem at each step in the recursion, and is (ii) guaranteed to generate a sequence of feasible dispatch policies for the RAC-OPF problem with nonincreasing costs. The proof of Proposition 3 can be found in Appendix B-C.

Proposition 3**.**

Let $v_{0}\in\mathbf{C}^{n}$ be a feasible solution to the zero-recourse AC-OPF problem (12), and define the matrix $V_{0}\coloneqq v_{0}e_{1}^{T}$. Let $\{x_{t},V_{t},\Lambda_{t}\}_{t=1}^{\infty}$ denote the sequence of solutions generated by the recursion in (13). The following properties hold for each step $t$ of the recursion.

1. (i)

Nonemptiness: $\mathcal{F}(V_{t})\neq\emptyset$. 2. (ii)

Cost montonicity: ${\rm tr}(MV_{t}^{H}A_{0}V_{t})\leq{\rm tr}(MV_{t-1}^{H}A_{0}V_{t-1})$.

We also remark that it is also possible to establish convergence of the recursive method (13) to a stationary point of the nonconvex problem ${\cal P}_{\rm II}$ using existing techniques from the literature [35, 36].

VI Case Study

We now illustrate the effectiveness of the proposed optimization method on the WSCC 9-bus [32] and IEEE 14-bus [34] test systems with different levels of renewable resource penetration and uncertainty. We refer the reader to the aforementioned references for the complete specification and single-line diagrams of the test systems considered. All modifications made to the original WSCC 9-bus and IEEE 14-bus test systems are summarized in Tables I and II, respectively. Additionally, for each test system considered, we permit load shedding at non-generator buses (i.e., load buses), where a reduction in load relative to the specified level is penalized according to the value of lost load (VOLL).666We refer the reader to Appendix D for a complete specification of the RAC-OPF problem with load-shedding. We set the VOLL equal to $4,000/MWh.

VI-A Renewable Generator Model

The real-time generating capacity of renewable generators $i\in\{2,\dots,6\}$ represents the only source of uncertainty in the power system being considered. Accordingly, we set the dimension of the random vector to $k=6$, and let the $i^{{\rm th}}$ element of the random vector $\bm{\xi}$ represent the maximum active power available to generator $i$ in real-time. In other words,

[TABLE]

for $i=2,\dots,6$. It will be convenient to our numerical analyses in the sequel to express the random vector $\bm{\xi}$ as an affine function of a zero-mean random vector $\bm{\delta}$ that is uniformly distributed over a unit ball. We define this relationship according to

[TABLE]

where the random vector $\bm{\delta}$ is assumed to have support

[TABLE]

It follows that the random vector $\bm{\xi}$ has support given by

[TABLE]

Here, $\mu\in\mathbf{R}^{k}$ and $\sigma\in\mathbf{R}_{+}$ represent location and scale parameters, respectively. In the following study, we set $\mu_{i}=15$ MW for each renewable generator $i\in\{2,\dots,6\}$. Qualitatively, the larger the scale parameter $\sigma$, the larger the a priori uncertainty in the real-time generating capacity of the renewable generators. The location and scale parameters are chosen in such a manner as to ensure that $\bm{\xi}_{i}$ respects the nameplate active power capacity limits for each renewable generator $i$ specified in Table I. We also require that $\mu_{1}=1$ to maintain consistency with our original uncertainty model in Section II-B. Finally, under the assumption that $\bm{\delta}$ has a uniform distribution, it is straightforward to show that the random vector $\bm{\xi}$ has a second-order moment matrix given by

[TABLE]

Renewable energy resources, like wind and solar, employ power electronic inverters, which can produce and absorb reactive power. The limits on the maximum and minimum amount of reactive power that can be injected by a renewable generator are determined by its inverter’s apparent power capacity, which we denote by $s_{i}^{\max}\in\mathbf{R}_{+}$ for each renewable generator $i$. It follows that the real-time complex power injection of each renewable generator $i$ must satisfy a capacity constraint of the form

[TABLE]

As the slight oversizing of a renewable generator’s apparent power rating is standard in practice, we set $s_{i}^{\max}=1.05\real\{g_{i}^{\max}\}$ for each renewable generator $i$. In order to ensure that Assumption 2 is satisfied, we enforce a more conservative form of the real-time apparent power capacity constraint (15) by setting the real-time reactive power limits for each renewable generator $i$ according to

[TABLE]

For simplicity, we also fix the nameplate reactive power limits of each renewable generator according to ${\rm Im}\{g^{\max}_{i}\}={\rm Im}\{\overline{g}_{i}(\bm{\xi})\}$ and ${\rm Im}\{g^{\min}_{i}\}={\rm Im}\{\underline{g}_{i}(\bm{\xi})\}.$ The reactive power limits specified in (16) specify the range of reactive power injections that are guaranteed to be available to a renewable generator in real-time, regardless of the active power supplied. Using the real-time active and reactive power capacity limits specified in (14) and (16), respectively, it is straightforward to construct matrices $\underline{G},\overline{G}\in\mathbf{C}^{n\times k}$ such that Assumption 2 is satisfied.

Finally, it is worth noting that the conservative generation limits $(\gamma_{i}^{\min},\gamma_{i}^{\max})$ defined in (11) admit closed-form expressions for the system being considered. They are given by

[TABLE]

VI-B Numerical Analyses and Discussion

We begin by examining the sensitivity of the generation cost incurred under the affine recourse policies that we propose to uncertainty in renewable supply. We do so by varying the scale parameter $\sigma$ from 0 to 15 in increments of 1.5, while keeping all other problem parameters fixed. It is worth noting that for $\sigma=0$, there is no a priori uncertainty in the renewable supply, and the RAC-OPF problem (8) reduces to the zero-recourse AC-OPF problem (12). In addition, for $\sigma=0$, we are able to verify that the solutions we compute are in fact an optimal solution by using the semidefinite relaxation of the AC-OPF problem described in [7]. For each value of $\sigma$ that we consider, we calculate an affine recourse policy according to the recursive algorithm specified in Eq. (13).777All numerical analyses were carried out in Matlab and semidefinite programs were solved using SDPT3 [37]. The machine used to to solve the problems has a 3.1GHz Intel dual-core with 16GB of RAM. We initialize the recursion with a feasible solution to the zero-recourse AC-OPF problem (12), which we compute using the Matpower interior point solver [33]. The recursive algorithm terminates when either one of the two following conditions hold: (i) the difference between the optimal values of successive iterations is less than $10^{-4}$, or (ii) the total number of iterations exceeds 500.

We plot the expected generation cost (and empirical confidence intervals) incurred by the affine dispatch policy returned by algorithm (13) versus the scale parameter $\sigma$ for the 9-bus and 14-bus systems in Figures 1(a) and 1(b), respectively. First, notice that, for each test system, the expected generation cost increases monotonically with the scale parameter. Such behavior is to be expected, as larger values of $\sigma$ correspond to larger uncertainty sets $\Xi$. It is also worth noting the ‘spread’ in the cost distribution induced by the dispatch policies that we compute also increases with $\sigma$. That is to say, renewable energy resources with a large variance in their real-time generating capacity will result in a larger variance in total generating costs. Such behavior is a consequence of the risk neutrality inherent to the expected cost criterion that we treat in this paper.

We examine the convergence behavior of the recursive algorithm (13) for the 9-bus and 14-bus test systems in Figures 1(c) and 1(d), respectively. For each system, we set $\sigma=7.5$ and plot (as dashed lines) the expected generation cost incurred at each step of the recursive algorithm (13) for five randomly generated initial conditions $V_{0}$.888We initialize the recursive algorithm (13) with randomly generated matrices $V_{0}$ by solving the zero-recourse AC-OPF problem (12) with with randomly generated marginal-cost parameters $\{\alpha_{i}\}$. Specifically, the generators’ marginal-cost parameters $\{\alpha_{i}\}$ are sampled independently and uniformly at random from the interval $[0,50]$. In both Figures 1(c) and 1(d), the solid black lines represent the expected cost trajectory returned by the recursive algorithm given the initial condition $V_{0}$ generated by the zero-recourse AC-OPF problem (12) under the true marginal-cost parameters specified in Tables I and II. It is interesting to note that the recursive algorithm converges to the same optimal value regardless of the initial condition. Additionally, the numerical results in Figures 1(c)-(d) are consistent with the guarantees of Proposition 3, which ensures that the recursive algorithm (13) will yield a sequence of feasible dispatch policies with nonincreasing costs. We also note that each iteration for the 9-bus (14-bus) system took 11.67 seconds (130.4 seconds) to complete on average.

In Figures 1(e) and 1(f), we plot empirical confidence intervals for the total active power load that is shed at non-generator buses under the affine dispatch policy returned by algorithm (13) versus the scale parameter $\sigma$ for the 9-bus and 14-bus systems, respectively. The confidence intervals are generated using 10,000 independent realizations of the random vector $\bm{\xi}$ for each value of $\sigma$. The lower and upper whiskers of the confidence intervals extend to the minimum and maximum values of the total active power load shed. For the 9-bus system, the affine dispatch policy never sheds active power load, as depicted in Figure 1(e). That is to say, the allowance of load shedding does not result in any reduction of load under the affine dispatch policies that we compute. For the 14-bus system, we empirically observe that only a small percentage of total active power load is shed under the affine dispatch policies that we compute. In particular, the total active power load shed never exceeds 0.16% of the total active power load at non-generator buses. Lastly, we remark that for both the 9-bus and 14-bus systems, we empirically observe that the affine dispatch policy never sheds reactive power load.

VII Conclusion

In this paper, we formulate the robust AC optimal power flow (RAC-OPF) problem as a two-stage robust optimization problem with recourse. Under the assumption that there is adjustable generation or load at every bus in the power transmission network, we provide a technique to construct a convex inner approximation of RAC-OPF in the form of a semidefinite program. In particular, the inner approximation is obtained by: (i) restricting the set of admissible recourse policies to be affine in the uncertain variables, (ii) approximating the semi-infinite constraint set by a sufficient set of finitely-many constraints, and (iii) approximating the indefinite quadratic constraints by majorizing convex quadratic constraints. Its solution yields an affine recourse policy that is guaranteed to be feasible for RAC-OPF. In addition, we provide an iterative optimization algorithm that generates a sequence of feasible affine recourse policies with nonincreasing costs.

There are several interesting directions for future research. First, affine recourse policies are likely to be suboptimal for the RAC-OPF problem. Thus, it would be interesting to investigate the design of convex relaxations for RAC-OPF to enable the tractable calculation of lower bounds on the optimal value of RAC-OPF. Such lower bounds can, in turn, be used to bound the suboptimality incurred by the feasible affine policies proposed in this paper. Second, the convex inner approximation technique developed in this paper relies explicitly on the rather limiting assumption that there is adjustable generation or load at every bus in the power transmission network. In order to accommodate the treatment of more general power systems, it will be important to relax this assumption, while preserving the robust feasibility guarantees developed in this paper. It would also be of interest to extend the techniques developed in this paper to accommodate discrete decision variables (e.g., unit commitment decisions) in the RAC-OPF problem.

Appendix A Concise Reformulation of RAC-OPF

Define matrices $\Phi_{i},\Psi_{i}\in\mathbf{H}^{n}$, for all $i\in{\cal V}$, and a matrix $E\in\mathbf{R}^{4n\times 2n}$ as follows:

[TABLE]

In addition, let $f\in\mathbf{R}^{4n}$ be a vector given by

[TABLE]

The RAC-OPF problem (8) can be reformulated as follows:

[TABLE]

Appendix B Proofs

B-A Proof of Lemma 2

Part (i):   Notice that

[TABLE]

The desired result follows from the fact that $A^{-}_{i}\preceq 0$.

Part (ii):   Let $N\in\mathbf{R}^{k\times k}$ be a Cholesky factor of $M$, i.e., $M=NN^{T}$. The existence of $N$ is guaranteed as $M$ is assumed to be positive definite. Since the matrix on the left-hand side of (19) is Hermitian negative semidefinite, we must have

[TABLE]

And since the trace of any Hermitian negative semidefinite matrix is nonpositive, we obtain

[TABLE]

The linearity of the trace operator and its invariance under cyclic permutations gives the desired trace inequality. ∎

B-B Proof of Proposition 2

Given a feasible solution $v_{0}$ to problem (12), we show that there exists $x_{0}\in\mathbf{R}^{2n}$ and $\Lambda_{0}\in\mathbf{R}^{m\times\ell}$ such that the point $(x_{0},V_{0},\Lambda_{0})$ is feasible for problem ${\cal P}_{\rm III}(V_{0})$, where $V_{0}=v_{0}e_{1}^{T}$. Let $x_{0}\in\mathbf{R}^{2n}$ be a vector given by

[TABLE]

where $S_{i}\in\mathbf{C}^{n\times n}$ is defined in (1). Since $v_{0}$ is a feasible solution to (12), it is easy to verify that $Ex_{0}\leq f$, where $E\in\mathbf{R}^{4n\times 2n}$ and $f\in\mathbf{R}^{4n}$ are defined in Appendix A. Therefore, it suffices to show that there exists $\Lambda_{0}\leq 0$ such that

[TABLE]

for all $i\in[m]$. Since $V_{0}=v_{0}e_{1}^{T}$, we have

[TABLE]

for all $i\in[m].$ Therefore, (20) takes the form

[TABLE]

for all $i\in[m]$. Fix an $i\in[m]$ and consider the constraint

[TABLE]

in problem ${\cal P}$. This constraint can be equivalently expressed as $v^{H}A_{i}v+b_{i}^{T}x\leq\pi_{i}$, where $\pi_{i}$ is the optimal value of the following convex quadratically constrained quadratic program

[TABLE]

In the above optimization problem, we have used the relation $\xi_{1}=1$ to rewrite $c_{i}^{T}\xi$ as $\xi^{T}C_{i}\xi$, where $C_{i}\in\mathbf{R}^{k\times k}$ is defined in problem ${\cal P}_{\rm II}$. By Asumption 1, there exists $\xi\in\Xi$ such that $\xi^{T}W_{j}\xi>0$ for all $j\in[\ell]$. This is a Slater condition, which guarantees strong duality to hold between the above optimization problem and its dual. Therefore, $\pi_{i}=\tau_{i}$, where $\tau_{i}$ is the optimal value of its dual problem, which is given by

[TABLE]

Since $v_{0}$ is a feasible solution to (12), it holds that

[TABLE]

Therefore, there exists $\rho_{i},\gamma_{i},[\Lambda_{0}]_{i1},\dots,[\Lambda_{0}]_{i\ell}\in\mathbf{R}$, such that $[\Lambda_{0}]_{ij}\leq 0,$ for all $j\in[\ell]$, $v_{0}^{H}A_{i}v_{0}+b_{i}^{T}x_{0}\leq\rho_{i}+\gamma_{i},$ and

[TABLE]

Since $v_{0}^{H}A_{i}v_{0}+b_{i}^{T}x_{0}\leq\rho_{i}+\gamma_{i},$ it follows readily that

[TABLE]

Subtracting $C_{i}+\sum_{j=1}^{\ell}[\Lambda_{0}]_{ij}W_{j}$ from both sides of (24), we observe that its left-hand side becomes equal to (22). Therefore, it suffices to show that

[TABLE]

as this implies that (22) holds. Using Lemma 3.1 in [38], the positive semidefiniteness constraint (23) can be equivalently described by the following set of conditions:

[TABLE]

Using the equivalent conditions in (26), we verify the satisfaction of the desired inequality in (25) in two separate cases: $\gamma_{i}=0$ and $\gamma_{i}<0$.

Case 1:  Let $\gamma_{i}=0$. Then, (23) implies that $\rho_{i}=0$. This follows from the fact that for any positive semidefinite matrix $X$, if the diagonal entry satisfies $[X]_{mm}=0$, then all off-diagonal entries in the correspoding row and column also satisfy $[X]_{ml}=[X]_{lm}=0,$ for all $l\in[n]$. In this case, (25) becomes equivalent to the first condition in (26), and we are done.

Case 2:  Let $\gamma_{i}<0$. By rearranging terms, the third condition in (26) holds if and only if

[TABLE]

It follows that (25) holds if

[TABLE]

For a matrix $X\in\mathbf{H}^{n}$, let $\eta(X)\in\mathbf{R}^{n}$ be the vector of eigenvalues of $X$ arranged in nonincreasing order. According to Corollary 7.7.4 in [39], it holds that $\eta(X)\geq\eta(Y)$ for any pair of matrices $X,Y\in\mathbf{H}^{n}$ that satisfy $X\succeq Y$. It follows that the third condition in (26) implies that

[TABLE]

where $\eta_{\max}$ is the largest eigenvalue of $C_{i}+\sum_{j=1}^{\ell}[\Lambda_{0}]_{ij}W_{j}.$ As this is a positive semidefinite matrix, we must also have that $\eta_{\max}$ is nonnegative. Therefore, $\rho_{i}\leq 2\sqrt{-\gamma_{i}\eta_{\max}}$. Using this inequality, we can bound the left hand-side of (27) from above. And therefore, it suffices to show that

[TABLE]

as this implies that (27) holds. Indeed, consider the function $f:\mathbf{R}_{+}\times\mathbf{R}_{++}\rightarrow\mathbf{R}$ given by $f(x,y)=-y+2\sqrt{yx}-x$, where $\mathbf{R}_{+}$ and $\mathbf{R}_{++}$ denote the sets of nonnegative and positive real numbers, respectively. To complete the proof, notice that the function $f$ has a maximum value of zero over its domain. ∎

B-C Proof of Proposition 3

For notational brevity, let us first define the matrices

[TABLE]

for each $i\in[m]$. In addition, define

[TABLE]

The proof is by induction on $t$.

Base step: Let $t=0$. By Propostion 2, there exists $x_{0}\in\mathbf{R}^{2n}$ and $\Lambda_{0}\in\mathbf{R}^{m\times\ell}$ such that $(x_{0},V_{0},\Lambda_{0})\in\mathcal{F}(V_{0})$. Hence, $\mathcal{F}(V_{0})\neq\emptyset$, and we established property (i) for the base step. Let

[TABLE]

It follows that

[TABLE]

since $(x_{0},V_{0},\Lambda_{0})\in\mathcal{F}(V_{0})$ is also feasible. It remains to show that ${\rm tr}(MV_{1}^{H}A_{0}V_{1})\leq J(V_{1},V_{0})$ in order to establish (ii) for the base step. Recall property (ii) of Lemma 2. Letting $V=V_{1}$ and $Z=V_{0}$, we obtain the desired inequality.

Induction step: Let $t=s$, and suppose that $\mathcal{F}(V_{j})\neq\emptyset$ for all $j<s$. We must show that $\mathcal{F}(V_{s})\neq\emptyset$ in order to establish (i). To do so, we use the assumption that $\mathcal{F}(V_{s-1})\neq\emptyset$, together with Lemma 2(i) to show that $(x_{s},V_{s},\Lambda_{s})\in\mathcal{F}(V_{s})$. Let $(x_{s},V_{s},\Lambda_{s})\in{\cal F}(V_{s-1})$ be an optimal solution to (13). This solution is guaranteed to exist since $\mathcal{F}(V_{s-1})\neq\emptyset$ by assumption. Thus, $Ex_{s}\leq f$, $\Lambda_{s}\leq 0$, and

[TABLE]

Fix an $i\in[m].$ By using Lemma 2(i) at $V=V_{s}$ and $Z=V_{s-1}$, we obtain

[TABLE]

Adding $-C_{i}+(b_{i}^{T}x_{s})e_{1}e_{1}^{T}-\sum_{j=1}^{\ell}[\Lambda_{s}]_{ij}W_{j}$ to both sides of (30), yields

[TABLE]

where the last relation follows from (29). Consider now the feasible set $\mathcal{F}(V_{s})$. We claim that $(x_{s},V_{s},\Lambda_{s})\in\mathcal{F}(V_{s}).$ Indeed $Ex_{s}\leq f,\ \Lambda_{s}\leq 0$ and $F_{i}(V_{s},x_{s},V_{s},\Lambda_{s})\preceq 0$ for all $i\in[m]$, where the last relation follows from (31) and the fact that $H_{i}(V_{s},V_{s})=V_{s}^{H}A_{i}V_{s}$. This completes the induction step for part (i).

We will now establish the induction step for part (ii). We must show that ${\rm tr}(MV_{s}^{H}A_{0}V_{s})\leq{\rm tr}(MV_{s-1}^{H}A_{0}V_{s-1})$ for any optimal solution $(x_{s},V_{s},\Lambda_{s})$ of (13). Let $(x_{s},V_{s},\Lambda_{s})$ be one such solution. Since it is optimal, we must have

[TABLE]

In the step of induction of part (i), we have shown that $(x_{s-1},V_{s-1},\Lambda_{s-1})\in\mathcal{F}(V_{s-1})$. Therefore, we obtain in particular

[TABLE]

It remains to show that ${\rm tr}(MV_{s}^{H}A_{0}V_{s})\leq J(V_{s},V_{s-1}).$ This follows by setting $V=V_{s}$ and $Z=V_{s-1}$ in Lemma 2(ii). ∎

Appendix C LMI Reformulation of ${\cal P}_{\rm III}(V_{0})$

For $i=0,1,\dots,m$, let $B_{i}=(A_{i}^{+})^{1/2}$ and $N=M^{1/2}$ be the square roots of $A_{i}^{+}$ and $M$, respectively. These matrices are guaranteed to exist since both $A_{i}^{+}$ and $M$ are positive semidefinite. In our reformulation, it will be convenient to write

[TABLE]

where $\text{vec}(\cdot)$ denotes the linear operator vectorizing matrices by stacking their columns. Let

[TABLE]

denote the part of $H_{i}(V,V_{0})$ that depends affinely on $V$. Applying the Schur complement formula to the matrix inequalities in ${\cal P}_{\rm III}(V_{0})$, we obtain the following equivalent reformulation of ${\cal P}_{\rm III}(V_{0})$ as a semidefinite program:

[TABLE]

Appendix D RAC-OPF with Load Shedding

In what follows, we provide an alternative formulation of the RAC-OPF problem, which allows for load shedding at non-generator buses in the power network. Here, a reduction in load relative to the nominal demand level is penalized according to a suitably chosen value of lost load (VOLL).999For example, the current VOLL used within the Midcontinent Independent System Operator (MISO) markets is $3,500/MWh [40].

In order to develop this generalization of the RAC-OPF problem, we first require some additional notation. Let $\mathcal{V}_{G}\subseteq{\cal V}$ denote the subset of buses connected to generators, and define $\mathcal{V}_{L}\coloneqq{\cal V}\setminus\mathcal{V}_{G}$ as the subset of non-generator buses. Also, let $n_{G}:=|\mathcal{V}_{G}|$ and $n_{L}:=|\mathcal{V}_{L}|$ denote the number of generator and non-generator buses, respectively. Define the complex load that is shed at each non-generator bus $i\in\mathcal{V}_{L}$ by $\lambda_{i}(\bm{\xi})$, where $\lambda_{i}\in\mathcal{L}^{2}_{k,1}$ is a recourse function determined by the ISO for each load $i\in\mathcal{V}_{L}$. The power balance equation at each non-generator bus $i\in\mathcal{V}_{L}$ can therefore be expressed as

[TABLE]

Finally, letting $\beta\in\mathbf{R}_{+}$ denote the VOLL, we arrive at the following reformulation of the RAC-OPF problem with load shedding:

[TABLE]

Appendix E Table of Widely-Used Variables and Symbols