key: cord-0128432-3f8j2hjw
authors: Liu, Haitian; Guo, Ye
title: Distributed Multi-Area Optimal Power Flow via Rotated Coordinate Descent Critical Region Exploration
date: 2021-10-29
journal: nan
DOI: nan
sha: 17fd7c098b5fba543dc5dc141a39b97fb917c923
doc_id: 128432
cord_uid: 3f8j2hjw

We propose a novel distributed algorithm, namely, the RCDCRE method. It originates from primal decomposition and tailors to the multi-area OPF. RCDCRE leverages the CRE method for privacy-preserve and fast convergence. Cyclic BCD is integrated to make RCDCRE more suitable for asynchronous updates in wild multi-areas. The coordinate system of RCDCRE is also rotated strategically to ensure convergence. We illustrate the simulation results through IEEE benchmarks.

Multiple utilities usually jointly operate a vast interconnected power system due to technical reasons and geographical partitions. As each utility is responsible for optimizing local assets, interchange scheduling is crucial for the overall efficiencies. For example, NYISO has an annual load level of 17.1 [GW] in 2020, where 13% of them come from the imports and exports substantial amounts of power from four adjacent control areas: New England, PJM, Ontario, and Quebec [1] .

Improving existing power interchange schemes contributed to clean energy innovation, providing solutions to climate change. With the growing trends of global warming, decarbonization is a pressing goal all over the world. In 2019, New York also passed legislation requiring 70% renewable in total power consumed by 2030 and 100% carbon-free by 2040. Increasing the total capacity of renewable energy resources (RESs) is the leading path towards decarbonization. However, the opposite geographical distributions of RESs and load centers require coordination among all the areas. For instance, the southwest of America has rich solar radiation, and midwest America has abundant wind sources according to NERL's data [2] . As mentioned, operators in northeast America like NYISO have significant net imports, yet they relatively lack RESs. Thus, scheduling inter-area dispatch can also reduce overall carbon emissions among these utilities.

Effective inter-area interchange offers potential abilities to reduce RESs' curtailment. In 2020 IEA's report [3] , CAISO curtails over 318 [GWh] RESs output in April 2020, reaching new record curtailment levels which are 67% more than in 2019. Two reasons cause the phenomena, 1) California had added 1[GW] of wind and solar PV capacity from May 2019 to March 2020, raising RESs output. 2) There is load-dropping of around 8% (mainly as a result of Covid-19-related measures) at the same time as solar output peaked. Hence, it is essential to schedule the interfaces to coordinate multiple utilities' flexibility. Such flexibilities can also provide resilience hedging H. Liu and Y. Guo RESs generation volatility and mutual support under extreme weather conditions. Ideally, the optimal utilization of tie-lines is the solution to a centralized OPF. Because multiple utilities jointly operate the interconnected power system. A distributed solution to the centralized OPF is preferred. The main challenge of distributed OPF is the trade-offs among convergence speeds, privacies, synchronization requirements, and coordination.

Active power is the primary concern in multi-area dispatch, which suits DC approximation. DC-OPF reforms the OPF into linear or quadratic programming, preferred in large-scale solving and distributed settings. Current methods tailored for power system operation mainly derived from dual decomposition, see [4] , [5] for a comprehensive review. The pros for such methods are convenient for implementation at the coordination level. The cons, however, are the convergence issues. As dual decomposition methods minimize local Lagrangian problems with relaxed constraints, all the intermediate solutions are not feasible until convergence. Moreover, an extra feasibility recovery step is generally necessary for dual decomposition [6] . Divergent cases may occur for ADMM methods [7] .

On the other side, primal decomposition has also made progress in power system applications. Typical methods includes marginal equivalent decomposition (MED) [8] , column and constraint generation (C&CG) [9] , and critical region exploration (CRE) [10] , [11] . The merit of such methods is the finite convergent property in certain problems. For example, MED identifies the marginal variables to exploit the finiteness of the structure of the active constraint set. C&CG is an extension of Benders decomposition. Rather than merely submitting cutting planes, such a method also exchanges binding variables to the coordinator to improve convergence. CRE leverages multi-parameter programming theory, which uses parametric mapping to preserve privacy. The geometric search of limited critical regions has ensured finite step convergence. There are no relaxed constraints when solving the subproblems with a fixed boundary state. Thus, empirically, such methods have better feasibility guarantees and convergence. However, privacies may be an issue as some methods have exchanged actual costs and network information (e.g., MED, C&CG). Although CRE is well privacy-preserved, it still suffers from the strong coordination with synchronization requirements. This paper proposes an alternative distributed solution to DC-OPF in terms of primal decomposition. Our methods leverage CRE for privacy-preserve and fast convergence. The synchronization requirements are removed by introducing an asynchronous block coordinate update scheme. Moreover, we also alleviate direct control of all coupling variables of a central coordinator, allowing each area to maintain its boundary phase angles independently.

II. MULTI-AREA OPF FORMULATION We use the DC power flow to approximate the network. For a multi-area system, the areas are coupled through tie-lines. The formulation is given by

where the variables include each area's power generation g i , internal phase angles θ i , and boundary phase angles θ i . The objective (1a) is to minimize the sum of each area's generation costs with coefficients Q i , c i . The network balance is subdivided into internal and boundary buses' power balance, shown in constraints (1b)-(1c). The notation j ∈ nbd(i) includes both area i itself and area j that are physically connected to area i. The line flows is reformulated into internal and tie-line flow constraints (1d)-(1e). Constraints (1f)-(1h) denote variable bounds. Note that there are no boundary generator for practical reasons. Constraint (1i) artificially assigns a reference bus, such that the solution to DC-OPF can generally be unique.

The selection of such reference bus is arbitrary.

To simplify notations, we eliminate (1i) and define

By reforming (1b), (1c), (1d), (1f), and (1g) into (2b), while (1e), and (1h) become instances of (2c). A standard primal decomposition model coupled by θ i , i = 1, . . . , N is given by

Before introducing the proposed method, we first review some preliminaries in multi-parametric programming and the critical region exploration (CRE) method. Each area's local problem parameterized in θ i , i = 1, . . . , N is given as

Critical regions (CRs) and value function (VF) are defined as 1

where A i,k is the active set index for area i at iteration k. The mapping relations are shown in [˜] values. For all active sets index that can optimally solve (3), we have

is piecewise affine/quadratic and convex. The feasible set Θ can be divided into a collection of CRs.

CRE aims to find a particular parameter that can minimize the sum of all local value functions. It comprises two steps 1) Find the optimal value for a specific VF and CR, given by

2) Do a sub-gradient search from current CR to adjacent one. Algorithm 1 summaries the implementation details.

Here, we use ∂J ⋆ , N Θ to indicate the sub-gradient set and normal cone w.r.t current optimal parameter θ ⋆ at iteration k. The following lemma follows from [10, Theorem 1].

Lemma 2. CRE method optimally solves (3) in finite iterations, and the θ ⋆ , J ⋆ opt are also optimal solutions of (2).

The limitation of the CRE method is that all parameters need to be managed by a central coordinator and updated synchronously. Such a method is generally hard to implement for a multi-area power system due to 2 practical reasons. 1) Coordinator has full control of all boundary states, 2) Clock synchronization requirement for multi-areas.

An intuitive idea to resolve such issues is to let each region optimize its boundary phase angles in sequence. Block coordinate descent (BCD) algorithm naturally suits the framework. We can conduct such iterations to minimize J ⋆ until reaching a stationary point. However, the dilemma is that the convergence of constrained BCD relies on the following two assumptions. Assumption 1. The objective J ⋆ is differentiable, if not, then the non-differentiable part is N block-separable [12] .

Assumption 2. The feasible space Θ is the direct sum of N separable real Hilbert spaces (H i ) 1≤i≤N [13] .

According to Lemma 1, J ⋆ is a piecewise convex function. Since feasible region Θ is usually a general convex polyhedron, these two assumptions cannot ensure. The BCD algorithm may be stuck at some non-optimal point. This paper proposes an alternative procedure, which we called rotated coordinate descent critical region exploration (RCDCRE). As the name indicates, the basic idea is to add an extra rotation operation when stucking. More specifically, we want to rotate the current parameter coordinate system into a particular position such that the convergence assumptions of constraints BCD hold. RCDCRE integrates CRE, BCD, and rotation. Algorithm 2 summarizes the implementation details.

if s is an odd number then t ← t + 1 ; Step 3-13, Algorithm 2 perform a single BCD on area i to minimize the value function J ⋆ . The sub-gradient search in CRE has been replaced by a coordinatewise descent, as indicated in Step 5. Note that the role of CRE here is not only to improve convergence, but also to preserve privacy. We remark that BCD cannot directly implement on the problem (2) without CRE, as variables x i and constraints (2b) are local private information that cannot be exchanged.

Step 15-19, we maintain a sub-gradient set in the outer loop. Multiple sub-gradients maybe added, depending on how many critical regions each area has explored at θ ⋆ . The BCD update is also realized in these steps, where we adopt a classic cyclic version, as indicated in Step 16 and 19. Readers could also try other asynchronous BCD variants in [12] .

Termination test for RCDCRE largely resembles to CRE, as shown in Step 21-26. However, notice that in RCDCRE the sub-gradient v is no longer used for searching but for realizing rotation. The coordinate system rotation proposed in this paper is denoted in Step 23 of Algorithm 2, which is briefed as generation of a specific rotation matrix R. To illustrate how to rotate the coordinate, let us first assume we have already obtained such a matrix. Then we can define the rotated parameter as

Consequently, the parameter coordinate system can be equivalently treated as a rotated one if we substitute θ = R −1 θ R into (2). The blocks in θ also retain the separation in θ R , namely, all area still control its boundary information. The rotation procedure completes by reform the problem (2), as indicated in Step 24 of Algorithm 2. Note rotation matrix R is orthogonal, the inverse operation is cheap since

Two things are of concern to conduct rotation, i.e., where to rotate the current coordinate system and how to realize such a process? Although Algorithm 2 can apply to any rotated coordinate system. There is no way to ensure convergence if we rotate arbitrarily because the new coordinate system may violate the BCD convergence assumptions. Therefore, in this paper, we proposed a strategic way to generate the rotation matrix R, which is summarized in the following proposition Proposition 1. Suppose current optimal parameter θ ⋆ is stucking (i.e., v ≥ ε optimal ). Then one of the 2 conditions must occur by rotating v in Step 21 of Algorithm 2 2 into any span of sub-block parameters θ i , i = 1, . . . , N 1) Finding a better parameter with J ⋆ (θ ⋆ new ) < J ⋆ (θ ⋆ old ), 2) Finding at least 1 new sub-gradient w.r.t current θ ⋆ .

The proof of this proposition originates from Lemma 2, since v in RCDCRE and CRE both indicate the same particular sub-gradient. By following the same arguments in [10] , the next exploration in direction v is guaranteed to find a new CR. Notice under Proposition 1, such direction has been rotated into the block coordinate descent direction in RCDCRE. Then we can assert that either there is a better parameter in the new CR or a new sub-gradient that can help us find another new CR after the next rotation, which completes the proof.

The last puzzle of the rotation procedure is how to generate a rotation matrix R that can rotate the parameter from θ to θ R . This procedure can be equivalently considered as to rotate between any high dimensional vector, which can be done in the following proposition.

Proposition 2. For any unit vector x 1 , x 2 ∈ R n , the rotation between x 1 and x 2 can be achieved by finite composition of Givens matrices.

Note that the Givens matrix generalizes 2D rotation into the following n dimensional rotation in the ij plane

where α is an appropriate rotation angle. By removing the i, j rows and columns, R ij becomes an identity matrix, which serves as rotation reference in the high dimensional 2D rotation. Obviously, any high dimensional rotation can be treated as a successive operation of 2D rotation. So with properly chosen Givens matrices, we can realize any rotation we want naturally. Under Proposition 1, we rotate sub-gradient v to a random vector in the span of sub-blocks. A simple choice of such random vectors is the standard basis. The specific rotation is realized by Proposition 2. Such a scheme can guarantee our next round exploration after rotation contributes to optima.

We illustrate how our proposed coordinate system rotation technique avoids the stucking of BCD. For ease of presentation, it is an oversimplified example that does not involve the CRE procedure. Interested readers can refer to [10] , [11] for illustrations on CRE. Our problem is given by minimize x1,x2

subject to

where there are two agents with an initial start point (−1, −1). It's easy to tell that the optimal solution to this simple problem is (−0.5, −0.5). As shown in Figure 1a , suppose x 2 do the coordinatewise minimization first, the optimal solution is (−1, 0). However, x 1 cannot find any improvement at (−1, 0), due to the problem violates the Assumption 2, namely, block separable structure of constraint sets. Notice that at the stucking point (−1, 0), a descent direction is induced by the set plus of normal cone and gradient. As the objective is differentiable, it is equivalent to the projected gradient given by

Note that [1, 0] T is unit anti-gradient of objective (6a) at (−1, 0). Here P is a projection matrix given by

where K is the collection of normal cone vectors. For this case, K is [1, 1] T , which is the same as the normal vector to constraints (6b). By normalizing the results in (7), we get the unit descent direction √ 2/2, − √ 2/2 T . Such procedure is a simplified version to the Step 21 of Algorithm 2. Under Proposition 1, we can rotate current descent direction into [1, 0] T , which corresponds to the standard basis of x 1 .

The transformation is given by

where the matrix in the equation is the rotation matrix R we wanted. As this is a simple 2D rotation, we don't need to refer to Proposition 2 to form such matrix. By substitute R into problem (6), the original problem is reformed as minimize x1,x2

as shown in Figure 1b . Current descent direction becomes [1, 0] T , the rotated problem can be solved by BCD readily.

In this toy example, we only show the rotation of stucking at a non-separable constraint set. Nevertheless, a similar analysis could also draw from stucking at non-separable and non-differentiable objectives, along with the CRE mapping procedure. We leave such discussions for interested readers.

We use MATLAB 2021a to perform the simulations, where all subproblems are solved via MOSEK version 9.2.17. Network data are obtained from MATPOWER 7.1 [15] . For illustration, the entire parameter space Θ has been divided into CRs using MPT 3.0 [16] in some cases. However, we need to point out that it is unnecessary when applying the proposed RCDCRE method. We adopt a cold start in all cases.

Consider the two-area power system shown in Figure 3 that stitched IEEE 14 and 30 systems together with 2 tie-lines as shown. All tie-line capacities are set as 10 MW, and internal line capacities are 100 MW. Linear cost coefficients c i of the generators were perturbed to c i := c i • (0.99 + 0.02ξ i ), for i = 1, . . . , N , where entries of ξ i are independent N (0, 1) (standard normal) variables. After eliminating bus 9 in area 1 as reference bus, the problem can be reformed into (3) with 2 agents and 2 dimensional parameters.

Simulation results that neglect quadratic generation costs are shown in Figure 2a . We can see that the proposed RCDCRE method can find the exact optimal solutions after switching the coordinates four times. The parameter space Θ contains 3 CRs, and we only need to explore two of them for the entire RCDCRE procedure. If we do not neglect the quadratic costs, then Θ comprises 4 CRs. We can see from Figure  2b , the problem stuck at a non-optimal point after switching coordinate two times. After the rotation, RCDCRE can find the optima by another two switchings, shown in Figure 2c . Note that the trajectory before the rotation is no longer parallel to the current axis.

We next show the results on a large four-area 472-bus network with 50 MW tie-line and 500 MW internal line capacities. The topology is depicted in Figure 4 . A quadratic generation cost is adopted, where the linear coefficients are perturbed similar to the two-area one. We compare the objective convergence of the proposed RCDCRE with three other approaches, i.e., CRE method taken from [10] , Consensus ADMM taken from [5, Algorithm 2], and Benders decomposition taken from [17, Section 5.1]. The results are shown in Figure 5 . CRE method and the proposed RCDCRE converge to the exact solution in 26 and 401 iterations, respectively. Meanwhile, it takes 712 and 589 iterations for Benders and ADMM to satisfy the relative optimality gap 3 . RCDCRE allows each area to update its boundary states rather than the coordinated update like CRE and Benders. Moreover, it enjoys a better convergence over ADMM due to adequate information exchanged. 3 We set the relative optimality gap as 10 −5 . 

This paper proposes an RCDCRE method that integrates CRE, BCD, and coordinate system rotation. The proposed RCDCRE leverages CRE for privacy-preserve and fast convergence, BCD to remove the synchronization requirements and alleviate central coordination. Moreover, we offer a strategic coordinate system rotation technique to reform the problem dynamically. By rotating the conditional steepest descent direction into any span of coordinate blocks, the proposed RCDCRE avoids the stucking even the original problem is violating BCD convergence assumptions. Proposed RCDCRE is tested on a 4-area 472-bus network, showing convergence superiority over ADMM and Benders decomposition methods.

2020 State of the Market Report for the New York ISO Markets

NREL's geospatial data science research

Renewables 2020 -Analysis and forecast to 2025

A Survey of Distributed Optimization and Control Algorithms for Electric Power Systems

Toward Distributed/Decentralized DC Optimal Power Flow Implementation in Future Electric Power Systems

Primal convergence from dual subgradient methods for convex optimization

The direct extension of ADMM for multi-block convex minimization problems is not necessarily convergent

A Marginal Equivalent Decomposition Method and Its Application to Multi-Area Optimal Power Flow Problems

Solving two-stage robust optimization problems using a column-and-constraint generation method

On Robust Tie-Line Scheduling in Multi-Area Power Systems

Coordinated Multi-Area Economic Dispatch via Critical Region Projection

Coordinate descent algorithms

Parallel random block-coordinate forward backward algorithm: A unified convergence analysis

Conditional subgradient optimization: Theory and applications

Matpower: Steady-state operations, planning, and analysis tools for power systems research and education

Multi-Parametric Toolbox 3.0

Introduction to Stochastic Programming, ser. Springer Series in Operations Research and Financial Engineering