1. Introduction. In many optimization problems, data is not completely known in advance. One of the main approaches to deal with this lack of information is robust optimization (RO). It assumes that the data lies in a predefined set of scenarios and the constraints have to be satisfied for any realization of the data in that set. RO does not require any knowledge about the distribution of the uncertain data and is useful when the distribution is hard to estimate and feasibility for a certain set of parameters is important. In reality, some optimization variables represent decisions that must be made before the actual realization of the uncertain data while other variables can be adjusted after the uncertain data becomes known. To account for such situations in RO, two-stage adjustable robust optimization (ARO) was introduced in [5]. Two-stage ARO problem includes two types of variables: the first-stage variables that are fixed and the second-stage variables that may change depending on the uncertainty realization. We refer to these variables as control and state variables, respectively, since this setting is typical for optimal control problems. ARO gives rise to more flexible decisions than robust optimization and thus could be less conservative.

Several approaches have been developed to solve ARO problems, most of them are approximations since solving ARO to optimality can be NP-hard even if the original problem without uncertainty is a linear program [5]. The hardness comes from defining the relations between the first-and second-stage decision variables, called decision rules, which are usually not specified and must be optimized. In certain problems occurring in practice, the exact functional form of the second-stage rules is predetermined but unknown since it is defined implicitly via, for instance, a system of non-linear equalities that are challenging, or even impossible, to eliminate. Implicitly defined decision rules are rarely considered in the literature while they occur in many applications, some of which are described below. This paper aims to partly close this gap and studies general quadratic ARO with decision rules implicitly defined by a system of quadratic equalities.

With rapidly increasing uncertainties in both the demand and supply of resources, networked infrastructure problems are important applications of ARO with implicitly defined decision rules. Examples of relevant networked infrastructures include electricity (see [8]), natural gas (see [1,30]), and water (see [42]). Operators must ensure that these systems remain in acceptable states despite uncertainties while also considering performance criteria such as operating costs. ARO provides a way to balance these potentially competing concerns, as discussed in [48] and [30]. As an illustrative application of ARO, we have chosen the so-called robust AC optimal power flow (ACOPF) problem that provides minimum cost operating points for electric power systems. ACOPF can be formulated as a non-convex quadratic optimization problem, as discussed in [19]. Even in the absence of uncertainties, ACOPF is NP-Hard (see [10]), and solving such problems under uncertainty for instances of realistic size is still a challenge, which we explain in detail in Section 5.

1.1. Contribution. This paper presents two main contributions. First, we propose a convergent iterative solution approach for two-stage non-convex quadratic problems under uncertainty. The resulting solutions are feasible for the underlying problem without uncertainty (in other words, the nominal problem) and approximately feasible for the ARO version of the problem. Second, we implement the proposed approach for ACOPF, and it outperforms the benchmark methods.

Solution framework for quadratic ARO under ellipsoidal uncertainty with implicit decision rules defined by quadratic equality constraints. The equalities increase the difficulty of ARO. Therefore, we address equality and inequality constraints separately. We construct piecewise affine approximations of the implicit decision rules. In particular, we express the state variables as functions of the control and uncertainty variables using the first-order Taylor approximations of the original implicit rules on small subsets. To our knowledge, we present the first algorithm for ARO which uses implicit second-stage rules. For each piece of the piecewise affine approximation, we eliminate the equalities and second-stage variables and obtain a standard non-linear quadratic problem in the first-stage variables under ellipsoidal uncertainty. This problem is reformulated into a semidefinite program (SDP) with (possibly) quadratic constraints. We suggest using an alternating projections algorithm to find locally optimal solutions for this SDP in the presence of quadratic constraints. Such a first-stage solution is feasible for the nominal problem, and possible constraint violations under uncertainty are limited by the parameters of our method.

Implementation for ACOPF and comparison with benchmarks. We apply the proposed framework to ACOPF with uncertainty in power supply and demand. The ACOPF is formulated as a non-convex quadratically constrained quadratic problem, and an ellipsoidal uncertainty set is considered. As shown in Section 5, our numerical results for ACOPF demonstrate the effectiveness of our approach in comparison to other approaches from the literature on small to moderate-size instances ranging from 6 to 118 buses, with the potential to consider larger instances.

1.2. Existing ARO solution approaches. One of the most popular ways to solve general ARO is an approximation where the second-stage decision variables can be written as affine functions of the uncertain parameters. This approach is proposed by [5], and its current state-of-the art is discussed in [17]. Besides affine decision rules, one could consider piecewise constant decision rules. Such rules are constructed by partitioning the uncertainty set into subsets and implementing a fixed second-stage decision for each subset. Naturally, the two types of decision rules could be combined into piecewise affine decision rules, which is suggested in [35]. Fixing the form of decision rules could restrict flexibility, so one could miss the optimal ARO solution. As an alternative, there exist convergent relaxations which gradually add constraint violating uncertain scenarios to the problem [9,6,49]. These approaches work best if the uncertainty set is polyhedral and the nominal problem can be solved efficiently, which is different from the setup in our paper.

Our approach has features of both described ARO frameworks. On the one hand, we use piecewise affine decision rules. However, those cannot be arbitrary good rules, as they must reflect an actual process represented via a system of equalities to prevent second-stage infeasibilities. Moreover, the equalities are not solvable for some first-stage decisions and uncertainty values. Thus, our method constructs piecewise rules that are close to the real ones when those exist and indicates a lack of second-stage solutions otherwise. Using approximations might still lead to second-stage infeasibilities, but one can control those by increasing the number of affine pieces. This behavior resembles the approaches ensuring feasibility on a subset of the uncertainty realizations.

Finally, let us discuss the literature addressing non-linear ARO with implicitly defined decision rules. Recently, some progress has been made in methods suitable for such ARO, and we outline the differences between our approach and these methods. First, we do not consider any assumptions on convexity and concavity, which makes our work different from all results in convex and linear ARO. Another advantage is that our approach, while using SDP, is not based on SDP approximations of the original problem, such as in [1,23,28,43,46]. We allow keeping some non-linearities from the original problem as opposed to classical SDP approaches that linearize all constraints. Additionally, independently of how much non-linearity one keeps, the resulting SDP constraints are of the size of the number of uncertain parameters. On the contrary, other approaches use SDPs of at least the size of the number of first-and second-stage variables. The number of uncertain parameters is often substantially smaller than the number of variables. For example, in ACOPF the number of variables depends on the number of buses in the system, uncertainty usually occurs in loads or renewable generators, and not all buses have those. Hence, the size of SDP constraints in our approach could be substantially smaller than in the above-mentioned papers.

Our approach differs from other robust optimization methods where non-linear constraints are linearized, such as [29]. We linearize equalities only, linearize locally, and our approximations are closely related to the original constraints via Taylor series. The results in this paper are also distinct from [32] since in the latter work a robust solution is obtained by iteratively tightening the inequality constraints. In [21], the authors tackle general non-linear optimization problems but use an alternative ARO formulation and thus a distinct solution strategy. Our approach is close in spirit to the approach in [38] which uses local Taylor expansions too. However, the authors of the latter paper work with chance constraints, use full linearizations and address specifically ACOPF. Similarly, the approach in [7] is tailored to linear complementarity problems.

This paper is organized as follows. In Section 2, we present the problem formulation and motivate our solution approach. In Section 3, we describe the proposed algorithms in detail. In Section 5, we evaluate our approach on ACOPF instances. Finally, in Section 6 we present conclusions and discuss directions for future work.

2. Problem formulation and general framework. We begin this section with the notation. We denote the range of matrix A by R(A). We denote the space of n×n symmetric matrices by S n and for A, B ∈ S n , the trace inner product of A and B is denoted by A, B := trace(AB). We use the notation [n] for the set {1, . . . , n}. For a vector V of length n, we denote the i th entry of V by V i . We say that two continuous maps f and g on a compact set A ⊂ R n are ε-close to each other if sup x∈A f (x) -g(x) ≤ ε for some given norm • . We say that two vectors (or elements of a vector space) a, b are ε-close to each other if a -b ≤ ε for some properly defined norm.

In this paper, we deal with quadratically constrained quadratic problems (QCQP). Our approach generalizes to problems of higher degree as we explain later, but we focus on QCQP for simplicity. Now, let n ζ , n y , n x , m eq , and m in be natural numbers and consider semialgebraic sets S y ∈ R ny , S x ∈ R nx defined by quadratic constraints. Consider quadratic mappings f : R ny → R, G :

Problem (2.1) occurs frequently in dynamic systems optimization, which inspired this paper. In such problems, an operator sets up the values of control variables y, and the state variables x are determined afterwards according to a system of equations defining the equilibrium. Some examples of dynamic systems optimization are energy problems (ACOPF as described in Section 5, optimal power dispatch presented by [11]), water problems (the valve placement problem by [18]) or gas problems (passive gas network feasibility problem studied by [1]). Moreover, problem (2.1) describes the more general class of bilevel optimization problems where after setting the values of y in the first stage, the second stage variables are chosen from the set of optimal solutions of the secondstage optimization problem. The KKT optimality conditions of the second-stage problem can be written as a system of equalities, where the final set of second-stage variables x consists of the original second-stage variables and Lagrange multipliers of the second-stage problem. As a result, one obtains the so-called complementarity formulation of the initial bilevel optimization problem, which has the form of problem (2.1) if the initial problem was linear. A typical example of a bilevel optimization problem is the Stackelberg competition in economics. Some examples of bilevel problems in engineering can be found in [37,3]. More information about linear complementarity problems and robust optimization approaches for them can be found in [7].

Remark 2.2. Problem (2.1) allows for binary variables and absolute values. One can write binary constraints on a variable a as quadratic equality constraints a = a 2 . One can write the constraint a = |b| as a 2 = b 2 , a ≥ 0. An example of problem (2.1) that is the main use case in this paper is the ACOPF problem. It can be written as follows using the general notation above. min y,x y P y + p y + p 0

where y are active powers and voltage magnitudes on PV buses, and x are voltages in rectangular form, see Section 5 for the full problem formulation. Another example is the valve setting problem in water distribution networks from [18] with the following formulation:

where y are pressure heads and valve placement indicators, and x are flow rates and absolute values of flow rates.

In practice, given the uncertain demand and supply that is encountered in such applications, one often aims at solving problem (2.1) under uncertainty, which results in the next formulation. and for any ζ ∈ Ω there exists x such that the following holds:

3) is a two-stage ARO problem with the uncertain parameter ζ. The first stage happens before the uncertainty realization. At this stage, one assigns values to the control variables y ∈ S y . The second stage happens after the uncertainty realization. At this stage, one has to choose the best feasible value of the state variables x for the given uncertainty realization. The goal is to select a value for y ∈ S y such that there would be a feasible solution x in the second stage for any uncertainty realization ζ ∈ Ω. Any solution to problem (2.3) is feasible for the underlying nominal problem and robust against potential uncertainty. We obtain the nominal problem by setting ζ = 0 in problem (2.3). Next, we summarize main assumptions used in this paper. This assumption is for simplicity. In fact, the problem only has to be quadratic in ζ, the requirements for other variables are milder. Our approach would be still applicable if the problem were polynomial of higher degree in y; see Section 4 for more details. The inequality constraints could be general polynomials in x as well. If the inequalities have higher degree in x, we can obtain degree-two polynomials using variable substitution and increasing the number of state variables and equality constraints. We emphasize that this procedure does not influence the final size of the problems we solve since this size only depends on the number of control and uncertainty variables, i.e., y and ζ. (d) Ω is an ellipsoidal uncertainty set of the form:

where Σ is negative semidefinite. Ω has a non-empty interior. We need assumptions on the degree of ζ and the shape of Ω to efficiently eliminate ζ from the problem using the S-lemma [47]; see Section 4.1. We have chosen an ellipsoidal uncertainty set as the base case since it is convenient for our approach; frequently appears in the literature as being less conservative than, for instance, box uncertainty; and has interpretations from both robust and chance-constrained perspectives; see, e.g., [20,13].

Assumption 2.5 (Assumptions without loss of generality). (a) S x is defined by inequalities.

The assumption is w.l.o.g. since if the definition of S x contains equalities, they could be moved to (2.3a) as a preprocessing step. If S x becomes unbounded afterwards, a large ball constraint for x can be added to the problem. (b) n x = m eq and there are no redundant equality constraints in (2.3a).

The assumption is w.l.o.g. If n x ≤ m eq and there are redundant constraints, they could be detected and eliminated as a preprocessing step. If n x > m eq , then some state variables are free and could be added to the pool of control variables. (c) Objective f is amenable for optimization (e.g., convex quadratic or linear).

The assumption is w.l.o.g. since if f is not convex, we introduce the epigraph control variable and add the epigraph constraint to S y .

Our goal is to approximate the original problem by a problem in control variables y only, since they represent the actual decisions to implement. The following idea motivates our approach: if we could analytically solve equalities (2.3a) for the second-stage variables x, we would express x as a function of y and ζ obtaining the second-stage decision rule p : R ny×n ζ → R nx . Then problem (2.3) would be equivalent to the following problem: Problem 2.6.

That is, if a known second-stage decision rule exists, then we could substitute it in the problem and eliminate the equalities and state variables to obtain a classical (not adjustable) robust optimization problem. Clearly, if the second-stage variables are determined from a system of non-linear equalities, there might be no unique analytical expression for the decision rule in the problem. However, such expressions exist on small subsets of S y × S ζ under known conditions according to the implicit function theorem (see, e.g., [41]). In essence, the approach we suggest restricts the optimization problem to such subsets and replaces the implicit decision rule by its first-order Taylor approximation on each subset. We present the approach in detail in the next sections.

3. Piecewise affine approximations of equality constraints. Our approach looks at piecewise affine approximations of the second-stage decision rules x = p(y, ζ) in problem (2.6), which are given implicitly via the system of equalities (2.3a), using Taylor series. In addition to its simplicity, the main advantages of the Taylor approximation are its good fit for the original function around the approximation point and the possibility to construct it for an implicit function. For the construction, we need some Jacobians of the system of equalities (2.3a). Denote L(y, ζ, x)

x be the Jacobians of L with respect to y, ζ, and x, respectively, evaluated at (ŷ, ζ, x). By Assumption 2.4(b), the Jacobian

x is a square matrix of the size n eq × n eq . Theorem 3.1 (The implicit function theorem: [41]

From here on we call the main condition mentioned in Theorem 3.1 the Implicit Decision Rules (IDR) condition for brevity: Definition 3.2 (The Implicit Decision Rules (IDR) condition). We say that the IDR condition holds for ŷ ∈ S y , ζ ∈ Ω if there exists x ∈ R nx such that L(ŷ, ζ, x) = 0 and the Jacobian J ŷ, ζ,x x is non-singular. We say that the IDR condition holds on sets Ŝy , Ω if it holds for each y ∈ Ŝy , ζ ∈ Ω.

We assume that it is possible to check quickly if the IDR condition holds for given ŷ, ζ since one can solve equalities (2.3a) rather efficiently in practice. This paper was inspired by engineering problems, such as ACOPF. For these problems, Newton's methods [50] are usually well-developed, fast, and precise enough to work for realistic problem instances. Theorem 3.1 implies that, if the IDR condition holds for some (ŷ, ζ), a unique local second-stage decision rule x = p(y, ζ) as in (2.6) exists on some open set around (ŷ, ζ), and the first-order Taylor approximation of this rule is

The above result applies to ŷ ∈ S y for which the Jacobian with respect to x is non-singular. This condition is in general desirable in applications to dynamical systems, which are the main target of this paper. For such applications, solutions with singular Jacobians may be physically unstable, see, for example, [45].

The idea is to work on small subsets of S y and Ω and use piecewise affine decision rules based on the above Taylor approximations within subsets. Fix some subset of S y and ŷ ∈ Ŝy . For M ζ > 0, partition Ω into Ω 1 , . . . , Ω M ζ and pick some ζ1 ∈ Ω 1 , . . . , ζM ζ ∈ Ω M ζ . If Ŝy is small enough and the IDR condition holds for ŷ and ζ1 , . . . , ζM ζ , the implicit decision rules from Theorem 3.1 exist on Ŝy and some open subset of Ω. Now, in problem (2.6), we impose the inequalities for each Ω k substituting the implicit rules with the corresponding approximation (3.2). The procedure results in an almost standard ARO approximation with piecewise affine decision rules (see problem (3.3) in Algorithm 1). The differences are that (i) y is restricted to a subset, and (ii) the decision rules are known, we do not need to optimize over them. We can solve the approximation and obtain some solution y * . As the last step, we could check that y * is feasible for some selected values of the uncertainty that are desirable in practice, for instance, that it is nominally feasible. Then we can repeat the procedure considering other subsets of S y and choose y * with the best objective value.

Intuitively, if the subsets of S y and Ω are small enough, the above approach should be able to discard solutions ŷ ∈ S y for which no decision rules from Theorem 3.1 exist and select solutions for which these rules exist. We develop this idea in the next corollary. k=1 Ω k by finitely many compact subsets such that the unique decision rules from Theorem 3.1 exist on products of Ŝy and each Ω 1 , . . . , Ω M ζ .

Proof. We begin by proving part (a). In each Ω k , k = 1, . . . , M ζ there is at least one ζ for which the IDR condition holds and therefore equalities (2.3a) have solutions. Hence Ω does not contain fully one of Ω k . To cover a ball of radius up to r with balls of radius in such a way that none of the smaller balls is fully inside the larger ball, we need to have r ≤ 2 . Thus we cannot have a violating subset that contains a ball of radius larger than 2 . Next we prove part (b). By 3.1. Approximation algorithm for general problems. We formalize the ideas from the previous section in Algorithm 1, and the approximation guarantee is presented in Theorem 3.4.

Algorithm 1: Piecewise affine approximations of problem (2.3) using Taylor series

where pŷ i , ζk ,x i,k (y, ζ) is defined in (3.2). Let y * ,i be the optimal solution to (3.3).

The next theorem shows that Algorithm 1 provides approximation guarantees and is able to find solutions robust in a "strong" sense (such that the inequality constraints hold with a margin).

Theorem 3.4 (Approximation guarantees of Algorithm 1 for problem (2.3)). Consider a δ > 0, and let y * = ∅ be a solution returned by Algorithm 1.

(a) There exists > 0 such that, if S i * and each subset in the partition Ω =

k=1 Ω k fit in balls of radius , then the following conditions hold on Ω, except for possibly a subset where no ball of size larger than 2 fits: (i) equality constraints (2.3a) hold, and (ii) each inequality constraint (2.3b), (2.3c) is violated by at most δ. (b) Let y * be not feasible for problem (2.3) on Ω ⊂ Ω whose Lebesgue measure constitutes at least β < 1 of that of Ω. If M ζ ≥ log θ log(1-β) and ζ1 , . . . , ζM ζ are sampled uniformly i.i.d. on Ω, the probability that y * passes the check in step 9 of Algorithm 1 is not larger than θ. Guarantees for feasible solutions (c) Consider the set S δ y ⊆ S y of all solutions y to problem (2.3) which (i) satisfy the IDR condition on Ω and (ii) are such that for any ζ ∈ Ω constraints (2.3b), (2.3c) hold strictly with a margin δ (i.e., the inequalities " ≥ δ" hold). Then there is a subset of S y and a finite cover of Ω for which Algorithm 1 will return an optimal y * ∈ S δ y . Proof. First, consider statement (a). We know that step 2 in the algorithm did not fail for y * and the corresponding partition Ω = M ζ k=1 Ω k . If ε is small enough, then by Corollary 3.3 (a) there are decision rules x = p(ŷ i * , ζ) on Ω, except for possibly a subset where no ball of size larger than 2 fits. A small enough ensures also that S i * y is such that the decision rules are valid for both ŷi * and y * . Thus part (i) of statement (a) follows. Now, consider inequality constraints (2.3b) (2.3c). By assumption (2.4) (c), these constraints are continuous, and there are finitely many of them. Since we work on compact sets, Taylor approximation of each decision rule pŷ i * , ζk

y and Ω k so that the approximated inequality constraints (3.3b) in problem (3.3) are δ-close to the original inequality constraints (2.3b), (2.3c) in problem (2.3). Thus the original constraints cannot be violated (within the subsets of Ω where the rules exist) by more than δ, so y * possesses property (ii). To prove statement (b), notice that with i.i.d. uniform sampling, the probability that all constraints of problem (2.3) are satisfied for one sample ζ is at most 1 -β. Hence the probability that they are satisfied for M ζ samples is at most (1 -β) M ζ , and if M ζ ≥ log θ log(1-β) , then the latter probability is not larger than θ. Finally, consider statement (c). By Corollary 3.3, there exist a compact subset of S y that contains y * and a compact finite cover of Ω where the implicit decision rules used by Algorithm 1 are valid. By the above argument for part (ii) of statement (a), if we make the subset of S y and each subset in the cover of Ω small enough, the Taylor approximations are so close to the actual decision rules that the inequality constraints cannot differ from their approximations by more than δ, and thus Algorithm 1 will choose y * if we restrict the search to subsets from S δ y . To our knowledge, Algorithm 1 is the first ARO algorithm which uses the second-stage policies defined by implicit functions, represented by equalities (2.3a). Theorem 3.4(b) speaks about uniform sampling from Ω, which is a numerically efficient procedure since Ω is an ellipsoid [14].

The main limitation of the algorithm in practice is the need to choose the subsets of S y and Ω. This is not trivial since inefficient partitioning may lead to long running times and numerically unstable problems. Constructing good subsets is a topic for separate research. For instance, to split Ω, one can use the results from [35]. In this paper, we are especially interested in testing Algorithm 1 as a proof-of-concept for robust ACOPF. This problem has an additional structure: control variables y appear linearly and separately in L. In the next subsection, we introduce a version of Algorithm 1 that uses this property and allows working on larger subsets of S y .

3.2. Approximation algorithm for problems where L is linear in the control variables. In this subsection, we assume that the control variables y appear linearly and separately in the equality constraints L(y, ζ, x) = 0 . That is, we can write

where L 1 is a polynomial function and J y is a constant matrix. As a result, the Jacobian of L(y, ζ, x) with respect to y is J y , and other Jacobians do not depend on y. Denote the Jacobians with respect to ζ and x for any y and some ( ζ, x) by

x , respectively. Then the first-order Taylor approximation

Observe that the Jacobian J ζ,x x and the above Taylor approximation are independent of ŷ. That is, the same approximation could be valid for various subsets of S y as long as the decision rules for these subsets and some Ω k ⊆ Ω lead to the same subset Ŝx ⊂ R nx . Hence it could be beneficial to use subsets of S x instead of subsets of S y , reconsidering Algorithm 1 as presented in Algorithm 2.

Algorithm 2: Piecewise affine approximations of problem (2.3) where L is linear in y (j) :

The IDR condition holds for some ŷ ∈ S y in ζk with solution xj ,

where pζ k ,x j (y, ζ) is defined in (3.5). Let y * be the optimal solution to problem (3.6).

The next Corollary 3.5 shows that the performance of Algorithm 2 is similar to the one of Algorithm 1. y , and look at the complement of this union. The distance between this complement and S j,k y is positive since they are disjoint, S j,k y is compact and the complement is closed. Hence, there is some > 0 such that if γ ≤ , then the restriction on the Taylor approximation in (3.6a) (**) pushes y * to be in the earlier constructed open cover of S j,k y , so there exist second-stage rules for y * around ζk . Moreover, we can adjust γ and and the size of Ω k to make sure that the Taylor approximations of the inequality constraints are δ-close to the original constraints, for any δ > 0. These are the conditions we used to prove Theorem 3.4 (a), thus the corresponding result holds for y * . Finally, statement (b) follows from the proof of statement (a) under the following consideration. If M ζ = 1, then statement (a) applies to the whole Ω = Ω 1 . Thus, if τ and γ are small enough and we obtain y * = ∅, the conditions (i) and (ii) of Theorem 3.4 (a) are satisfied on the whole Ω.

Algorithm 2 could be more efficient than Algorithm 1 since it could capture more values of y ∈ S y , even if it looks around only one solution x. However, we cannot guarantee that a strictly feasible solution y * as defined in Theorem 3.4 (c) will be found under some conditions as we cannot control y as precisely as in Algorithm 1. Also, to say that at least one of the conditions (j) in (3.6a) holds for each Ω k , one may have to introduce additional binary control variables.

By Corollary 3.5 (b), if Ω is small enough for the chosen approximation error, it suffices to find one good small subset of S x , and there is no need to partition Ω. The nominal problem (2.1) for which we want to find a robust solution is usually feasible and, moreover, its optimal solution is usually "stable" in the sense that the Jacobian in (3.6a) is non-singular. Thus, for small uncertainty sets, one could start from some nominal feasible control solution ŷ1 at hand (e.g., the original optimal solution), find the corresponding state variables solution x1 such that L(ŷ 1 , 0, x1 ) = 0 and restrict Algorithm 2 to work around x1 as in (**) in (3.6a). If the algorithm finds a solution ŷ2 = ŷ1 , one could set x2 such that L(ŷ 2 , 0, x2 ) = 0 and repeat the procedure trying to find a better solution. Notice that a similar procedure would be valid for Algorithm 1 too, by Theorem 3.4 (a), but we would construct subsets around y 1 , y 2 in S y instead of subsets around x 1 , x 2 in S x . Now, we have constructed two algorithms, one is suitable for general problems (2.3) (Algorithm 1), and one is tailored to problems (2.3) linear in control variables (Algorithm 2). To solve the optimization problems (3.3) and (3.6) in the algorithms, we need to deal with quadratic inequality constraints under ellipsoidal uncertainty, which is the topic of the next section. The rest of this section is needed if the initial problem (2.3) contains inequalities that are nonlinear in y or ζ. If it turns out that problem (3.3) in Algorithm 1 or problem (3.6) in Algorithm 2 is linear in y and ζ, then classical robust optimization techniques apply to the corresponding problem. In this case one can use a standard reformulation of a linear constraint under the ellipsoidal uncertainty (see,e.g., [5]) and rewrite the problem as a second-order cone program. Finally, from here on, we focus on problem (3.6) for brevity, but the method is valid for problem (3.3) as well.

4.1. Eliminating the uncertainty from the problem. Under Assumption 2.4, all inequalities in problem (2.3) are of degree at most two in x and ζ, hence so is h in (4.1):

for some given parameters A, B, b, c, d, and a function g which contains all monomials that are non-linear in y. The precise values of the parameters and the form of g depend on the functions in (3.6a) and can be obtained directly from those functions. Recall that Ω has the form (2.1). Such a combination of h and Ω allows reformulating (4.1) and eliminating the uncertainty from it due to the S-lemma. Then the following two statements are equivalent:

S-lemma is well known in robust optimization (see [4]) but is usually applied to convex problems in y, ζ while we use it for a general quadratic constraint (4.1).

Proposition 4.2. Constraint (4.1) with h as in (4.2) and the uncertainty set (2.1) holds if and only if there exist λ, γ such that

where Σ, σ and r are the parameters from (2.1), and other parameters are defined in (4.2).

Proof. We have

If some inequalities in problem (2.3) are linear in x and ζ, then h in (4.1) is linear in ζ as well. In this case, the constraint considered in this section simplifies to the classical linear robust constraint. It is well known how to eliminate uncertainty from such constraints; see, e.g., [4]. Each linear constraint under the ellipsoidal uncertainty in (2.1) will become a second-order cone constraint after eliminating the uncertainty. By eliminating ζ, we reformulate problem (3.6). This reformulation can be written in a general way as Problem 4.3.

where γ are auxiliary control variables as in Proposition 4.2, y includes other original and auxiliary control variables, H is a polynomial mapping, F is a linear mapping, and C is a proper semialgebraic cone amenable for optimization (e.g., positive semidefinite cone, second-order cone, or the nonnegative orthant).

Under Assumption 2.4, the equality constraints are quadratic, the objective function is convex with degree at most two, and C is the positive semideninite cone. The approach we use next is applicable for general polynomial equality constraints and semialgebraic cones C, therefore we introduce the formulation above.

Problem (4.3) would be a classical conic problem for which many solvers exist if we did not have those non-linear equality constraints. Hence, we first relax those complicating equalities and obtain a solution feasible for (4.3c). Then we iteratively transform it into a solution feasible for the whole problem (4.3). Notice that problem (4.3) cannot be unbounded as the feasible set of its nominal problem is compact by Assumption 2.4. We begin with a lower bound for the problem, which can be obtained using any relaxation of the polynomial equality constraints. One can use lifting techniques, e.g., SDP relaxations or the reformulation-linearization technique as in [40]. The lower-bound relaxation might provide a solution that is infeasible for (4.3b). To obtain a feasible solution from the solution to the relaxation, we use the alternating projection method as presented in Algorithm 3. Let z U be an upper bound on f (y) in problem (4.3). Define two sets: We denote y := (y nc , y c ), where y nc , y c are the subsets of variables y that are involved (resp. not involved) in the non-convex constraints (4.3b). We use two subsets of the variables to speed up the algorithm by skipping the iterations in which only y c changes. Proof. Item (a) follows by construction of the algorithm since we limit the number of iterations. Item (b) follows from Assumption 2.4 and Theorem 7.3 in [15] as the sets A and B are closed and semialgebraic. Item (c) follows from the fact that A is convex. Hence the normal cone at any point in the interior of A equals {0}, and thus A and B satisfy the conditions of Theorem 2.1 in [15], which implies linear convergence. Next we prove item (d). For a given ν, let S ν := A ∩ B ∩ {(y, γ) : f (y) ≤ νf (y * )+(1-ν)z L }. Let (y i , γ i ) be the best solution of the algorithm at iteration i, and assume that it is not the closest local optimum denoted by (ŷ, γ). Then f (y i ) < f (ŷ) and there exist ν i > 0 such that S νi = ∅ and (y i , γ i ) is close enough to S νi to satisfy the conditions of item (b). Notice that the latter argument might fail if we require (ŷ, γ) to be a global optimum. Using ν i , the algorithm finds the next point (y k+1 , γ k+1 ) that is tol-close to S νi and such that f (y i ) < f (y i+1 ) ≤ f (ŷ). Since the interval (f (y i ), f (ŷ)) is bounded, for any δ > 0 there are N large enough, tol small enough and a sequence (ν i ) N i=1 such that f (y N ) is δ-close to f (ŷ). In the next subsection, we demonstrate the performance of the combination of Algorithm 2 with Algorithm 3 to solve problem (2.3) on the ACOPF. The experiments show that in the majority of tested cases this combination of algorithms finds a robustly feasible solution within 15 minutes of computational time for small uncertainty sets.

Optimal power flow (OPF) is one of the key optimization problems relevant to the operation of electric power systems. OPF solutions provide minimum cost operating points that satisfy both equality constraints termed the "power flow equations" which model the power system network and inequality constraints that impose limits on line flows, generator outputs, voltage magnitudes, etc. Accurately modeling the steady-state behavior of power systems requires the non-linear AC power flow equations, which can be formulated as a system of quadratic equality constraints.

Compounding the difficulties posed by the power flow non-linearities, rapidly increasing quantities of wind and solar generation are introducing significant amounts of power injection uncertainties into electric grids. To address these uncertainties, researchers have studied a wide range of stochastic and robust OPF problems [39], many of which use the DC power flow approximation; see, e.g., [44,8] for several relevant examples. This linear power flow representation permits the application of stochastic and robust optimization techniques developed for linear programs. Alternative approaches replace the AC power flow equations with other more sophisticated approximations, Adjust (y 1 , γ 1 ) := 1 tol (y 0 , γ 0 ) (to proceed with the while loop) Return "The best obtained solution is (y * , γ * )" such as the work in [34] and [38], or convex relaxations, such as the work in [43]. Such approaches can provide useful solutions in many contexts, particularly when the approximations are iteratively updated or adaptively adjusted. However, the quality guarantees from these approaches are provided with respect to the approximation or convex relaxation as opposed to the original non-convex ACOPF problem. Since power flow approximations and relaxations can introduce substantial errors relative to the non-linear AC power flow equations [36,16,2], the resulting solutions may lead to unacceptable constraint violations during operation in the physical system.

The power systems literature also includes approaches that directly address the non-linear AC power flow equations. These approaches can provide high-quality solutions in certain instances but may be limited to special classes of problems, such as systems that satisfy restrictive requirements on the power injections at each bus as in [29]. Other approaches use scenario-based techniques that enforce feasibility for selected uncertainty realizations, possibly obtained via subproblems that compute worst-case uncertainty realizations with local solvers as in [12] or convex relaxations as in [28]. Certifying robustness with such approaches is challenging due to the possibilities of local solutions and inexact relaxations. Rather than seeking the worst-case uncertainty realizations, the approach in [32] instead bounds the worst-case impacts of the uncertainties with respect to each constrained quantity. While this approach provides guarantees regarding the satisfaction on the engineering inequality constraints, each iteration requires the solution of many computationally expensive subproblems. We also note recent work in [26] that uses so-called "convex restriction" techniques (see [25]) to compute robustly feasible ACOPF solutions. While promising, this approach is undergoing continuing development and requires specialization to the particular non-linearities in each class of problems.

Accordingly, new computational methods such as the framework proposed in this paper are needed to address power system optimization problems that model both uncertain power injections from renewable generators and non-linearities from the AC power flow equations. The proposed framework would most naturally fit into a system operator's generation scheduling processes as part of near-real-time (five-minute to hourly) generator setpoint computations via optimal power flow algorithms, which is a key application at the heart of power system operations. Since optimization in this case is repeated frequently during the day, the uncertain deviations in loads or generation are likely to be of small size. However, they can lead to infeasibility of the nominal solutions, thus an approach is needed that takes the uncertainty into account. This situation corresponds well with the result in Corollary 3.5 (b) about small uncertainty sets. Extensions to incorporate binary variables would facilitate other applications, such as solving unit commitment problems for day-ahead scheduling that models generators' start up and shut down decisions. Thus, developing algorithms suitable for such extensions is an important direction for future work.

Exploiting the polynomial representation of the AC power flow equations, we next apply the approach described in this paper to the robust ACOPF problem, beginning with our notation and the problem formulation. Consider a power network with the set of buses N = {1, . . . , n} and the set of lines connecting these buses E. We denote the set of buses with generators by G and the active and reactive power demand (load) at each bus k ∈ N by P d k and Q d k , respectively and denote the index of the reference bus by s. To implement thermal restrictions on the transmission lines, we impose line current limits, see [50]. Our objective is to minimize the cost of power generation, which is one of classical objectives in OPF problems.

Denote the active and reactive power injections due to load or generation fluctuation by P r k and Q r k , respectively, for all k ∈ N . In the nominal ACOPF without uncertainty, P r k and Q r k are known and fixed. Next we define the ACOPF problem as a quadratic optimization problem. Problem 5.1.

where the last constraint sets the phase angle of the reference bus to zero. Now, let the active power fluctuations for each k ∈ N be P r k = P r k + ζ k , where ζ represents the uncertainty. We assume that the power injection uncertainties from the load and generation at each bus k ∈ N are modeled via a constant power factor cos φ k so that the reactive power fluctuations are

Without loss of generality, we let P r k = Qr k = 0, otherwise one can adjust the loads P d k and Q d k . We denote by δ the total change in the active power losses due to the redistribution of power flows from the uncertain power injection fluctuations relative to the losses from the nominal power flows. Note that δ is typically near zero, as the losses themselves are usually small and the changes in losses are even smaller. For an operating point to be robustly feasible, the generators must account for the total change in the active power injections, n i=1 ζ i -δ, associated with each uncertainty realization without leading to constraint violations. We adopt a "participation factor" model where this change in power injections is distributed among all generators according to a linear recourse policy with specified participation factors α k for each generator k. Thus, for each k ∈ N , the actual active power generation consists of the nominal power P g k and an adjustment in generation due to the uncertainty:

Thus, when introducing uncertainty to problem (5.1), we replace P g k in this problem with (5.1). We note that this model represents the steady-state behavior of widely used automatic generation control (AGC) (see [22]) and is adopted in many robust and stochastic OPF formulations, e.g., those used by [43], [38], and [32]. To model the uncertainty, we let the uncertain parameters ζ belong to the region Ω = {ζ ∈ R n ζ : ζ Σζ ≤ 1}, where Σ is a covariance matrix. That is, our uncertainty region is an ellipsoid centered on the point with no fluctuations.

For k ∈ N , we denote by V g k := x 2 k + x 2 k+n the squared voltage magnitude at bus k. Following traditional power system modeling practices, we consider three types of buses: PV, PQ and the reference bus. If k is a PV bus, the active power P g k and squared voltage magnitudes V g k are set by the operator while the reactive power Q g k can change. If k is a PQ bus, then the active power and reactive power are fixed to constant values while the voltage magnitude can change. Without loss of generality, we assume that active and reactive power generation at PQ buses is zero, otherwise the loads can be adjusted. Finally, the operator selects the voltage magnitude at the reference bus while the active and reactive powers are free to vary. We also introduce a variable t that denotes the worst-case upper bound on the active power on the reference bus. We use this bound to estimate the worst-case objective value over the uncertain power injection fluctuations, as is typical in robust optimization problems. As a result, the control variables y in the problem include t, P g k , where k belongs to the set of PV buses, and V g k , where k belongs to the union of PV and reference buses. Now we define the problem in the same form as (2.3) to more easily use the results from the earlier sections. This yields the following:

In the next subsection, we run numerical experiments solving problem (2.3) with the inputs defined above and the instances from Matpower [50].

In this section, we implement Algorithm 2 and compare it with several algorithms used in the literature. All computations are done using MATLAB R2021a and Yalmip (see [27]) on a computer with the processor Intel ® Core ® i7-8665U CPU @ 1.90GHz and 16 GB of RAM. Semidefinite programs are solved with MOSEK, Version 9.3.21 [33].

In the experiments, we need to choose the vector α and the parameters defining the ellipsoid Ω. We assume that all uncertainty balances out on the reference bus since this is the default setup of Matpower. Therefore α s = 1 and α k = 0 for k = s. We implement Algorithm 2 with the following inputs:

• We do not split the uncertainty set Ω and set ζ = 0.

• We consider only one x: the initial feasible solution to the solution of the nominal problem (5.1) provided by Matpower [50]. 3), linearizing all non-convex quadratic terms in the control variables and obtaining an SDP constraint of the size n y + 1. • In the ACOPF problem, some matrices A in (4.2) are negative semidefinite. We do not have to project on the corresponding equality constraints (4.3b). Instead, we replace the equality sign by "≤" and add the resulting constraints to the definition of A in (4.3). • We do the posterior check in step 2 of Algorithm 2 for ζ = 0, i.e., we check that the obtained solution is nominally feasible. Setting ν i = 1 implies that we do not try to improve the objective value after finding the first feasible solution. We considered setting ν i < 1, however, for all tested cases the initially obtained objective value was close to the lower bound. Attempts to improve the objective values resulted in using substantially more or all N iterations with negligible or no improvement in the objective value. This result can be explained by the high quality of the lower bound from the SDP relaxation, which frequently provides a feasible or close to feasible solution to problem (4.3). Therefore, we decided to stop at the first obtained feasible solution by setting

As benchmarks, we use three other approaches. First, we consider the nominal solution found by Matpower to see how robust or not it is. Next, we implement the classical DCOPF relaxation of ACOPF, which is a linear problem. We follow the approach from [8]. In this paper, the parameters α are optimized while in our case they are fixed. We can eliminate all equalities and second-stage variables in DCOPF, thus obtaining a standard linear robust optimization problem under ellipsoidal uncertainty. Finally, we implement an approach based on the classical SDP relaxation of ACOPF [24]. Without the uncertainty, this relaxation usually provides exceptionally good approximations. A number of papers in the literature (e.g., [43,46]) implement an ARO version of this relaxation by replacing the nominal positive semidefinite matrix in the relaxation by a linear matrix function of the uncertainty; see the earlier mentioned papers for details. Since the initial relaxation is linear, except for the SDP constraint, after substituting the linear decision rule, all equalities can be eliminated and the standard linear robust approach applies. To ensure approximate feasibility of the uncertain SDP constraint, one usually imposes this constraint for many enough realizations of the uncertainty. We followed the approach in [43] and imposed the SDP constraint for the largest possible value of each uncertainty. We have also tried using the sampling approach to ensure feasibility of the constraint with high probability, as in [46]. However, even for small cases this required many scenarios, resulted in too many SDP constraints, and lead to numerically unstable problems even if we used row generation to add the constraints. For all algorithms, we limit the running time to 15 minutes (for our approach, 15 minutes from the moment after obtaining the SDP lower bound). DCOPF and Taylor approach do not reach this time limit in the presented cases. For the SDP relaxation, this limit is reached for the large cases (57 and 118 buses), therefore this approach is omitted in the corresponding table.

The two benchmark models we consider have affine constraints and use a chance constraints framework. However, to obtain the final reformulations, both models exploit the analogy between the ellipsoidal uncertainty and chance constraints for normally distributed random variables and affine constraints. Therefore, we also use this analogy to conduct experiments. Namely, we say that the uncertainty vector ζ is a random variable with the mean zero and the vector of standard deviations is equal to w% of the initial load. Let ζ be normally distributed, and let W be a variancecovariance matrix of ζ. Let P be such that W = P P (obtained, e.g., via Cholesky decomposition) and q 1-ε be the (1 -ε)-percentile of the standard normal distribution. Then making one affine constraint robust against the uncertainty set Ω := {ζ : P -1 ζ ≤ q 1-ε } (5.7) is equivalent to requiring this constraint to hold with probability at least 1-ε% for ζ. This is a well-known result which follows from classical robust optimization (see, e.g., [4]) applied to a linear constraint with ellipsoidal uncertainty and a straightforward reformulation of a chance constraint using the cumulative probability function.

We allow the uncertainty to occur at all buses with positive active power loads. This is easy to change. For instance, we could also consider uncertainty in generation by allowing additional sources of uncertainty at generator buses. We use two options for the variance-covariance matrix W . For the first option, we set W = Diag (wP d k ), where Diag is the operator that creates a diagonal matrix from a vector. We let w vary from 0.01 (1% of the load at the corresponding bus) to 0.5 (50% of the load at the corresponding bus). In the second option, we allow random correlations among the uncertainties. Namely, we generate a random positive definite matrix and scale it such that the standard deviations are equal to wP d k . We relate our uncertainty to 95% chance constraints using Σ = P -1 P -1 , ρ = 0 and r = -1.65 2 in (2.1). In real-life applications, one can choose a correlation matrix that is suitable for the given application; for instance, one can assume that correlations are proportional to distances between buses.

Choosing the probability of violations of the chance constraints is not trivial. The number of constraints in the problem is large even for small cases. Thus, when we restrict the probability that each chance constraint is violated, the probability that at least one constraint is violated tends to one when the number of constraints grows. Formally, to obtain the 5%-safe radius of the ellipsoid for m constraints, we can just require all chance constraints to hold with probability 1 -ε/m. However, in practice this approach is usually too conservative, and thus we choose the typical 5% level for each constraint and investigate how robust the solutions are for this setup.

Another reason to connect robust constraints to a probabilistic setup is the property that the volume of the ball with a constant radius tends to zero when the dimension of the ball grows. Thus, for larger cases, it would be hard to generate a large enough sample of uncertainties within the unit ball while we can easily generate uncertainties from a normal distribution.

To evaluate the performance of the algorithms, we generate ζ according to the normal distribution around zero with the covariance matrix W . In particular, we first generate a standard normal vector z and then use ζ = P z, where P comes from the earlier mentioned factorization W = P P . We generate 1000 values of the uncertainly including the nominal case where ζ = 0.

The computational results are presented in Tables 123. Tables 1 and2 cover power networks with fewer than 30 buses, and Table 3 shows the results for larger power networks. Our main indicator is the share of experiments where the solution provided by the corresponding approach had some infeasibilities. We round all infeasibilities down to 1 × 10 -3 per unit since this is precise enough for our application. We also checked more precise rounding, and the differences in the tables were minor. In all presented instances, the equality constraints were always feasible.

For each case, we present all levels of uncertainty where at least one model has at most 10% of constraint violations. As we mentioned earlier, it is not obvious which performance should be considered good, especially for large cases. A perfect performance would be to have less than 5% of experiments with constraints violations since this is the violation possible in our setup for a single constraint. However, given that even the smallest case we present has more than ten constraints, we choose 10% as an acceptable violation.

We do not present typical Matpower smallest cases with three buses (LMBM3) and five buses (WB5) since for the first case even the nominal solution is robustly feasible for moderate values of the uncertainty, and for the second case none of the models is able to obtain a robustly feasible solution. The latter test case is known to be especially challenging in various contexts (see [31]).

In all tables, the first column denotes the values of the uncertainty as a fraction of load, mentioned earlier as w, in percent. The second column indicates the evaluated model. The first two model names speak for themselves, the third name "SDP" means the model from [43], and the name "Taylor" indicates the approach from this paper. The third column shows average estimated objective values among all instances. In general, robust optimization approaches do not aim to optimize the average objective, so the second column is just a rough indicator of the performance of each approach in terms of the objective value. We find this indicator sufficient since our primary goal is not to compare the objectives but rather to ensure feasibility of the constraints. The fourth column shows the running time of each benchmark algorithm required to find a solution used for the experiments. We do not include the input construction time for each test case since all benchmarks use roughly the same inputs. The fifth column highlighted by grey shows our main indicatorthe share of all experiments where the solution provided by the corresponding approach had some infeasibilities. The last two columns show the maximum number of violated constraints among all experiments: "PQ" indicates violations of active and reactive power bounds in one test case and "VI" indicates violations of voltage and current flow bounds in one test case. Small violations of the "PQ" constraints could lead to more serious problems than small violations of the "VI" constraints, therefore we have separated these indicators.

A hyphen indicates that the corresponding model delivered no solution (the resulting problem these approximations could become looser under uncertainty. Adding some tightening constraints to the "SDP" model could potentially improve the robustness of that approach.

As to the running times, our approach the slowest for small test cases since it is an iterative approach due to Algorithm 3. when approach found solutions, could usually do that within several of the alternating projections 3. Longer running times for the two largest instances indicate that the alternating projections algorithm went through many iterations and could not converge, thus a robust solution may not exist for the considered subset of state variables. For the two largest cases, the SDP approach is the slowest, and it is in fact too large to solve (e.g., for case 118 it includes 100 SDP constraints of the size 236 × 236).

Finally, our approach often provides somewhat higher average objective than others. Given that the approach is more robust than the others, the difference in the average objective is minimal. To understand the performance with respect to the objective better, we looked at Case 9 with correlations, which is robustly feasible for three of four models for the 1%-uncertainty and robustly feasible for two models for the 10%-uncertainty. Below we show the box plots of the objective values of all robustly feasible approaches. The two box plots use same interval of possible objective values in the y-axes for comparability. We see that the "Taylor" approach tends to have higher cost in both cases, and it also results in larger standard deviations for the larger uncertainty, meaning that very low costs are also more common than in DCOPF. The worst-case objective value of our algorithm in the experiments is a little larger than for DCOPF, and the algorithm does not find the DCOPF solution. To find more solutions, such as the DCOPF solution when it is robust, it could be useful to enlarge the search space of the algorithm.

Finally, we see that adding correlations to the uncertainty ellipsoid does not substantially change the results for small cases. However, when the instance size grows, no robust solutions could be found by any model for "case30", "case57", and "case118", even for the lowest values of the uncertainty. Therefore, we do not provide results with correlations for larger cases. Further investigation is needed to see why exactly this situation occurs and to possibly tune the models to work for larger cases with correlated uncertainties.

Our evaluations are equally conservative for all models. First, we say that more than 10% of experiments with constraint violations is a bad performance while we have many constraints and the analogy between the normal distribution and our ellipsoid is only valid for one constraint. The 5%-safe radius of the ellipsoid for many constraints would be much larger, but using this radius would be too conservative. Second, we do not deeply analyze magnitudes of violations. By the experiment setup they are larger than 1 × 10 -3 per unit, but the importance of such violations can depend on the type of constraint and the details of the test case. Finally, we did not fine-tune the algorithm to each specific case and allowed for as many power injection uncertainties as possible. The performance could improve if the algorithms are fine-tuned, especially for large instances.

6. Conclusions and directions for future research. In this paper, we propose a framework to obtain approximately feasible solutions to quadratic ARO problems with equalities, which implicitly define the second-stage decision rules for the state variables. We replace the implicit decision rules by their explicit piecewise affine approximations. As a result, we can eliminate the state variables from the problem and replace the original ARO problem by a sequence of classical quadratic problems with additional tractable conic constraints. Since, a generally efficient algorithm does not exist yet for the latter problems, we design an alternating projections algorithm that converges to a local optimum of the problem. For any ε > 0, if the piecewise affine approximations are fine enough, the suggested algorithm guarantees that the second-stage equality constraints are satisfied and the inequality constraints are not violated by more than ε on "large" subsets of uncertainty, where "large" is defined in Theorem 3.4 (a), (b). The feasibility of the second-stage equalities is rarely addressed in the literature, thus we consider analysing and ensuring such feasibility an important contribution to the existing research. We suggest two versions of the algorithm, Algorithm 1 for general problems and Algorithm 2 for problems linear in the control variables.

We implement the algorithm for ACOPF problems with uncertainty in loads and simulate the uncertainty to evaluate the performance of the algorithm in comparison to three known benchmarks: nominal solution, robust DCOPF and robust SDP relaxation. The solutions provided by our approach are robustly feasible for small uncertainty sets and for cases with up to 118 buses. Moreover, the solutions are substantially more robust than the benchmarks. The algorithm also performs well in terms of the objective function values. The experiments show that DCOPF and SDP approximations work well for problems without uncertainty, but possible inexactness of these methods is magnified after adding the uncertainty. As a result, they are less robust in an ex-post experiment. In contrast, our approach yields more robust results in the experiments. The good performance can be explained by the fact that our approach preserves more non-linearities from the original problem. Notice that the size of the final problem (i.e., problem (3.6)) which we solve depends on the number of sources of uncertainty. Hence, even for larger networks, the approach could work well for systems with uncertain power injections at a limited subset of the buses.

We conclude by analyzing the limitations of our algorithms and suggesting directions for further research. First, the non-linearities which we keep in the problem have drawbacks. Namely, when DCOPF can find a robust solution, it does so substantially faster than the approach we propose. At the last step of our approach, after removing the state and uncertainty variables, we need to solve a problem with quadratic and SDP constraints. We proposed the alternating projections in Algorithm 3 for this purpose. This algorithm only finds locally optimal solutions, may take a number of iterations to converge, and each iteration solves an SDP problem, which makes the procedure rather slow for larger instances. However, this approach is still relevant since no generally efficient algorithms exist for such problems yet, and our Algorithm 3 contributes to the development of quadratic optimization with SDP constraints. In the future, one could use an alternative algorithm or a relaxation to solve the obtained problem with quadratic and SDP constraints.

Next, our algorithms can exploit various subsets of uncertainty, control and/or state variables. However, the actual implementation considered small uncertainty sets without partitioning them, and we restricted the search to one subset of state variables around the initial optimal solution. As expected, for larger values of the uncertainty in our numerical experiments, the algorithm becomes imprecise or could not find a feasible solution. A natural remedy to increase precision and find additional feasible solutions would be to consider several subsets of the state variables as described in subproblem (3.6). Next, to work with larger uncertainty sets, we could combine Algorithm 1 and Algorithm 2 with the approach proposed in [35] to split larger uncertainty sets.

Third, some improvements could be done in particular for ACOPF problems. These problems possess much sparsity, and it is to a large extent preserved under the transformation (3.2). Therefore, it would be also important to explore the true scalability of our method for ACOPF problems by implementing the algorithms more efficiently and investigating possibilities to exploit the inherent sparsity structure of ACOPF problems. Finally, it would be interesting to look at the extensions of our method for multiperiod problems, especially those incorporating binary variables to enable application to unit commitment problems for longer-term generator scheduling.