An Efficient Homotopy Method for Solving the Post-Contingency Optimal Power Flow to Global Optimality

Optimal power flow (OPF) is a fundamental problem in power systems analysis for determining the steady-state operating point of a power network that minimizes the generation cost. In anticipation of component failures, such as transmission line or generator outages, it is also important to find optimal corrective actions for the power flow distribution over the network. The problem of finding these post-contingency solutions to the OPF problem is challenging due to the nonconvexity of the power flow equations and the large number of contingency cases in practice. In this paper, we introduce a homotopy method to solve for the post-contingency actions, which involves a series of intermediate optimization problems that gradually transform the original OPF problem into each contingency-OPF problem. We show that given a global solution to the original OPF problem, a global solution to the contingency problem can be obtained using this homotopy method, under some conditions. With simulations on Polish and other European networks, we demonstrate that the effectiveness of the proposed homotopy method is dependent on the choice of the homotopy path and that homotopy yields an improved solution in many cases. For at least 5% of the test cases, bad local minima were identified, and the homotopy method yielded a solution that was significantly better than state-of-the-art interior point methods in terms of reducing the violation cost during a catastrophic contingency scenario.


I. INTRODUCTION
Optimal power flow (OPF) is a fundamental tool for power system network analysis, where the goal is to find a low-cost production of the committed generating units while satisfying the technical constraints of the system. Finding a global optimum for a large-scale OPF problem modeled with AC power flow equations is a difficult task due to its innate nonconvexity, but it is nevertheless crucial for the reliable and efficient operation of power systems.
Initiated by the work [2], conic optimization has been extensively studied in recent years and proven to be a powerful technique for solving OPF to global or near-global optimality. The paper [2] has indeed shown that a semidefinite programming (SDP) relaxation is able to find a global minimum of OPF for a large class of practical systems, and [3] has discovered that the success of this method is related to the underlying physics of power systems. [4] and [5] have developed different sufficient conditions under which the SDP relaxation of OPF provides zero duality gap. Moreover, [6] has found an upper This work was supported by grants from ARO, ONR, AFOSR, and NSF. Parts of this paper have appeared in the conference paper [1]. Compared with [1], the new additions to this paper are major theoretical results, analysis of generator outages, and comprehensive case studies on large-scale networks.
bound on the the rank of the minimum-rank solution of the SDP relaxation of the OPF problem, which is leveraged in [7] to find a near globally optimal solution of OPF via a penalized SDP technique in the case where the SDP relaxation fails to work. These ideas have been refined in many papers to improve the relaxations via branch-and-cut approaches, conic hierarchies, and valid inequalities [8]- [11]. The reader is referred to the survey paper [12] for more details.
Power operators are concerned with security-constrained OPF (SCOPF) instead of an idealistic OPF problem. SCOPF can be regarded as a large number of OPF problems coupled to each other via physical constraints, where the first OPF corresponds to the operating point of the system under the normal condition and the remaining ones are associated with contingencies. It is shown in [7] that SDP relaxations are able to obtain high-quality solutions of SCOPF. However, since SCOPF is a gigantic problem with an enormous number of variables, even simple local search methods are ineffective for real-world systems. The common practice is to approximate SCOPF as a single OPF for the base-case subject to many surrogate contingency constraints expressed in terms of the variables of the base-case. SDP and many other methods may be used to solve such problems, but they will only find the voltages for the base-case. It is essential for the operators to determine the values of the controllable parameters for each contingency in advance. This leads to a high number of decoupled optimization problems. The objective of this paper is to solve each of these problems, named contingency-OPF, to global optimality using fast algorithms.

A. Homotopy for Optimization
Homotopy methods have been used to improve the convergence of optimization problems. While convergence to a global minimum with probability one is guaranteed for a convex problem [13], this is generally not true for nonconvex problems. In order to understand when homotopy can be effective in finding a global solution for nonconvex optimization, we explore a minimization problem of the form: min x f (x) where f : R n → R is a nonconvex function of x ∈ R n . This problem is named (P o ). Note that the function f (·) can incorporate exact/inexact penalty functions to enforce constraints on x, implying that this formulation is general for both unconstrained and constrained optimization [14]. We will call (P o ) the "base-case" problem. A deformed version of the base-case, which is also a nonconvex minimization problem, will be denoted by (P f ) and defined as min xf (x). We consider two possible methods for solving the deformed problem that are based on local search: Fig. 1. Evaluating the performance of homotopy on one-dimensional unconstrained minimization problems. The figure compares two different samples (1) and (2), with two different methods (a) and (b). The dotted lines show how the solution from the previous iteration is used in local search algorithms to solve the next problem. The red dots show the solution at each iteration using the position of the dotted lines as the initial point. For the one-shot method (a), the solution of P o is used as the initial point for P f . For the homotopy method (b), the base problem P o is gradually transformed to P f over three iterations, updating the initial point as the solution to the previous problem. a) One-shot method: Use the solution of P o as the initial point for any descent numerical algorithm to solve P f . b) Homotopy method: Generate a (discretized) homotopy map from P o to P f . Use the solution of P o as the initial point, but update it at each step of the homotopy by solving an intermediate problem using local search that is initialized at the solution of the previous step. A linear (non-discretized) homotopy map can be defined as: with the property that P (0) = P o and P (1) = P f . Depending on f (x) andf (x), homotopy may or may not lead to better results than solving the deformed problem in one shot. In Figure 1, we see an example where homotopy is effective in finding the global minimum of a deformed problem and another example where it leads to a non-global local minimum. Knowing when homotopy will be effective is highly dependent on understanding how the shape of the function changes from the base-case to the deformed problem. In the current literature, there is a lack of theoretical results to characterize the performance of homotopy in finding a global optimum.

B. Contributions
In this paper, we develop a homotopy method to improve the quality of the contingency-OPF solution. Instead of solving for the solution to a contingency-OPF problem directly via a descent numerical algorithm, we generate and solve (using local search) a series of intermediate optimization problems wherein we gradually remove a component of the power system. We show that the effectiveness of homotopy to find a global solution of the contingency-OPF problem is dependent on the homotopy path, and therefore, we characterize desirable homotopy paths. In doing so, we prove that the contingency-OPF generically has a unique global minimum. Note that it is essential to find a global solution because constraint violations in the case of a contingency are very expensive to deal with and a global solution corresponds to the minimum violation.
In Section II, we present the formulation of base-OPF and contingency-OPF. Next, in Section III, we introduce the homotopy method that connects contingency-OPF to base-OPF via parametrization. In Section IV, we develop theoretical results to characterize cases when homotopy will lead to a global solution of the deformed problem. Finally, in Section V we implement the homotopy method on actual test cases and verify its effectiveness. The proofs and additional simulation results appear in the Appendix.

C. Notations
The symbol R N denotes the space of N -dimensional real vectors and (·) T denotes the transpose of a matrix. Re{·} and Im{·} denote the real and imaginary parts of a given scalar or matrix. The symbol | · | is the absolute value operator if the argument is a scalar, vector, or matrix; otherwise, it is the cardinality of a measurable set. The symbol denotes the elementwise multiplication between two vectors. Let 1 n and 0 n denote the n-dimensional vectors of ones and zeros, respectively. Furthermore, 1 k n denotes an n-dimensional vector of ones except for the k-th element that is zero. The imaginary unit is denoted by j = √ −1. Let the power network be defined by a graph G(V, E), where V is the node set and E is the edge set. We assume that there is one generator at each node. Each node i ∈ V has an associated complex voltage v i , a fixed demand P d i + jQ d i , an unknown generation p g i + jq g i , and we assume that the nodal shunt admittance is zero. The complex voltage v i can be expressed in polar form, v i = |v i |e jθi , where |v i | and θ i denote the voltage magnitude and phase angle at bus i, respectively. With a slight abuse of notation, |v| denotes the vector of all voltage magnitudes. In addition, we define θ ij = θ i − θ j . The set of neighboring nodes of node i is denoted by N (i). Each line connecting two nodes i and j is represented by a standard Π-model with a series admittance y ij = g ij + jb ij and a shunt admittance y sh ij = g sh ij + jb sh ij . Then, the nodal admittance matrix Y is defined as Finally, p ij and q ij are the real and reactive power flows from bus i to j, respectively.

II. FORMULATION OF THE DISPATCH PROBLEMS
In this section, we present the mathematical formulations of the security-constrained base-OPF and the contingency-OPF. The "classic" optimal power flow problem without security considerations is a static optimization problem: where f o (·) represents the operating cost, and the decision variables x = (|v|, θ, p g , q g ) ∈ R 4|V| represent the vector of voltage magnitudes, voltage phase angles, real power generations and reactive power generations. We assume that there is a generator at each bus because setting P min i = P max i = 0 would account for the case when there is no generator. The constraints model technical limits, such as the power flow equations and explicit bounds on variables. Nonlinearities are introduced into the constraints via AC power flow equations, and these nonlinearities along with the voltage magnitude lower bounds result in the nonconvexity of the problem. The classic-OPF problem can be expressed in a standard optimization form as follows (note that h(·) and g(·) are vectors):

A. Security-constrained Optimal Power Flow
Suppose that there is a set of possible contingencies, namely K, where each contingency corresponds to a line or generator outage. Each contingency k ∈ K introduces a new set of variables x k , and therefore, for a network with |V| buses and |K| contingencies, the SCOPF problem will involve optimizing over 4|V|(|K| + 1) scalar variables. The contingencies also add operational constraints of their own. In addition, there are physical limitations on how the post-contingency network can adapt from the base-case, and these limits are added as constraints that are functions of the base-case variables.
Since this extremely high-dimensional problem is cumbersome to solve, in practice the contingency constraints are approximated via methods such as LODF and PTDF [15]. In essence, this approximates the contingency variable x k as a function of the base-case variable x. Therefore, postcontingency equality constraints for contingency k are approximated by a composite function of the form h k (x) t k (a k (x)), where a k (x) represents the control actions that are taken in the event of a contingency. The same goes for postcontingency inequality constraints, represented by g k (x).
Finally, another important consideration is how SCOPF performs when the problem is infeasible. Therefore, we model some operational limits using soft constraints with extra variables that capture the amount of violation. The objective function that is minimized is the sum of real power generation costs in the base-case as well as a weighted sum of equality constraint violation penalties in the base-case and contingencies. The standard optimization form is presented below: where φ(·) and φ k (·) represent the penalty functions for the violations. We denote this SCOPF problem as the base-OPF, distinguishing it from the contingency-OPF presented next.

B. Contingency-OPF
Recall that the SCOPF solves for the base-case operating point by taking into account the possible failures in the network. In the process, it approximates the relationship between the contingency operation point x k and the basecase operating point x. However, it does not actually solve for the x k 's. Therefore, for each contingency we need to solve a contingency-OPF problem to find the best operating point for the specific contingency scenario, given the base solution. This problem resembles the classic-OPF problem, except that there are additional coupling constraints that tie the problem to the original base-case. For instance, the voltage magnitude at a bus should be equal to its base-case value unless the reactive capacity of the generators at that bus is exhausted.
We model a contingency, such as a line or generator outage, by changing the system parameters from their base values. For example, a line outage physically means that power cannot flow over that connection, which can be modeled by setting the resistance of the line to infinity (this changes the admittance matrix). In the event of a line outage, the power is re-routed through other paths and therefore the amount of loss in the system changes. However, the difference in loss is small enough such that there is no need for additional participation from other generators, unlike in the scenario of a generator outage. Therefore, we can fix the real power generation to be equal to the base-case values and solve for the remaining variables such that the violations for the bus balance equations are small and spread out as much as possible. This is because a large concentrated violation in a few buses can result in serious issues for the power network, whereas small power mismatches can be taken care of by real-time feedback controllers. Taking these into consideration, each contingency-OPF under study is given as where we have the pair of equality constraints above for every i ∈ V and the set Ψ ⊂ R 5|V| is defined as below: The set V q is the set of buses that hit their upper or lower reactive power generation bounds in the base-case, and |v i | base is the voltage magnitude of bus i in the base-case. The notations G ij , B ij reflect the potential change in the admittance. Note that real power generation is now a fixed parameter obtained from a solution of the base-OPF and therefore has been denoted by capital P g .
For generator outage contingencies, there is an additional aspect to consider. A generator outage corresponds to setting the real power generation at that generator to zero. However, in order to compensate for the lost generation, the system operator needs to increase the power generation at other generators that participate in the outage response. The above framework is general enough to incorporate this difference: where P g,f is the new setpoint for the real power generation. Denoting x = [|v|, θ, q g , σ p , σ q ] as the combined variable, contingency-OPF in a standard optimization form would be: Note that f (·) is not the same objective function used in classic-OPF or the base-OPF but merely a simplified notation for φ(σ p , σ q ). In this paper, we focus on the case when where c p i and c p i are cost coefficients. Similarly, h(·) is the not the same as the constraint functions used in classic-OPF or base-OPF.
If the optimal objective value of the contingency-OPF is zero, it means that the solution of the base-case could be modified to stay feasible in case of the contingency. However, the focus of this paper is on hard instances with a nonzero optimal cost, meaning that some of the constraints must be violated to accommodate the outage. In these cases, since taking corrective actions to deal with nodal power violations is expensive, it is essential to find a global solution.

III. METHODS
Homotopy and continuation methods have long been used in mathematics and engineering to solve systems of nonlinear algebraic equations [16]. Continuation methods in mathematics describe the continuous transformation of an easy problem into the given hard problem [17].
The benefit of homotopy methods compared to other iterative methods is that homotopy methods may yield global rather than local convergence. These methods are most useful for problems where convergence to a global solution is heavily dependent on a good initial point, which can be hard to obtain. More recently, probability-one homotopy methods have been applied to solving optimization problems. The applications include optimal control ( [18], [19]) and statistical learning [20]. Typically, the homotopy methods in optimization focus on parametrizing the KKT conditions ( [17], [21]) or the objective function ( [13], [22]). Our method is similar to the homotopy optimization method described in [13], wherein a series of local minimization problems are solved. However, we will focus on a more generalized theoretical analysis of homotopy, allowing for a homotopy map on the set of constraints.
Homotopy methods have also been applied in the field of power systems, primarily to solve the power flow (PF) problem for cases that do not converge. The continuation power flow (CPF) problem is used to find a set of solutions of the power flow problem, starting at some base load and ending at an operating point near the voltage stability limit [23]. Homotopy methods have been shown to improve convergence of the PF problem ( [24]- [26]) as well as compute all possible solutions to the PF problem [27]. In the following subsections, we present a homotopy method that parametrizes the contingency-OPF to model a gradual line or a generator outage.

A. Homotopy Method for a Line Outage
In order to solve the contingency-OPF problem, we propose a homotopy method that gradually changes certain parameters of the problem from the base-OPF, rather than abruptly changing the structure of the network. For a line outage contingency, we introduce an aggregate homotopy parameter λ = [γ, β, γ sh , β sh ] corresponding to the series admittance and the shunt admittance, where γ, β, γ sh , β sh ∈ R |E| . To be more precise, we parametrize the admittance in the contingency-OPF as follows: which creates a family of OPF problems, named H λ , which can be written in the standard form of: Now, let ∈ E be a line that connects buses i and j, and consider a contingency scenario in which the line is out. Notice that λ o = [1 |E| , 1 |E| , 1 |E| , 1 |E| ] corresponds to the original network before the line outage, and λ f = [1 |E| , 1 |E| , 1 |E| , 1 |E| ] corresponds to the post-contingency network after the line outage. By varying λ from λ o to λ f , the homotopy map allows us to create fictitious power networks that constitute a series of intermediate OPF problems. Hereby, let ρ(λ) = 0 represent a continuous homotopy path in the space of λ that starts from the point λ o and ends at the point λ f .

B. Homotopy Method for a Generator Outage
For a generator outage, our proposed homotopy map gradually decreases the real power generation at the generators that are out and gradually increases the real power generation at the generators participating in the contingency response.
For the simplicity of presentation, let us consider contingencies associated with a single generator (generator k) outage. This is common practice in power systems and is referred to as the N − 1 criterion. Yet, note that the proposed method can easily be extended to multiple generator outages and is incorporated in Algorithm 2.

Algorithm 1 Homotopy-OPF for Line Outages
Given: 1. Power network G(V, E) 2. Contingency set K with line outages l k ∈ E for k ∈ K 3. Discretized homotopy path Λ := {Λ 1 , . . . , Λ T } Initialize: Solve base-OPF problem to find a globally optimal solution (v 0 , σ, σ k ) and obtain p g * , σ p * , σ q * through calculation. Formulate the contingency-OPF problem: 1. Fix real power generation to base-case solution: P g = p g * 2. Find V q , the set of buses that hit their upper or lower reactive power generation bounds in the base-case. for k ∈ K do Define (i,j) as the from and to buses of line l k .
Let P g,o ∈ R |V| be the real power generated at all generators in the base-case. Using the participation factors of generators that are still active in the contingency, we can compute P g,f ∈ R |V| , the real power generated at all generators after the contingency. Since generator k is down in this contingency scenario, P g,f k = 0. A more detailed illustration of how to determine P g,f using participation factors will be available in the Appendix. Similar to what we did for line outage contingencies, we introduce an aggregate homotopy parameter λ = [γ, β] with γ, β ∈ R |V| to create the following homotopy map: Focusing on the first equation where we parametrize the real power generation, notice that λ o = [1 |V| , 1 |V| ] corresponds to the original network before the generator outage, and λ f = [0 |V| , 0 |V| ] corresponds to the post-contingency network after the generator outage. By varying λ from λ o to λ f , the homotopy map allows us to trace a gradual generator outage. Equation (5) parametrizes the reactive power demand, and we will set the value Q d,f Q d,o . Although the justification for this extra parametrization is not clear for the moment, we will explain later that the parametrization needs to be of high enough dimension in order for the homotopy method to be effective. The series of homotopy problems have the same form as those for the line outage, given by Equation (3).

C. Implementation of Homotopy-OPF
Starting with a solution to the base-OPF, we aim to iteratively solve a series of homotopy-OPF problems along the homotopy path to eventually arrive at the contingency-OPF. So far, we have described the homotopy path in the continuous sense. However, in practical implementations, we Algorithm 2 Homotopy-OPF for Generator Outages Given: 1. Power network G(V, E) 2. Contingency set K with generator outages R k ⊂ V for k ∈ K 3. Discretized homotopy path Λ := {Λ 1 , . . . , Λ T } Initialize: Solve base-OPF problem to find a globally optimal solution (v 0 , σ, σ k ) and obtain p g Formulate the contingency-OPF problem: Define P g r as the fixed real power generation at r ∈ V Define ∆P g k as the total lost real power generation at k: Find V q , the set of buses that hit their upper or lower reactive power generation bounds in the base-case. 2. Remove real power generation for generators in R k : P g r ← 0 ∀r ∈ R k 3. Compute participation factors α g r for r ∈ V \ R k 4. Add real power generation for participating generators: and P g,f r as the initial and final real power generation values at generator r for all r ∈ V Let P g,o ← p g * and P g, Solve homotopy-OPF H Λ i using initial point (ṽ,σ p ,σ q ) and obtain new solution (v, σ p , σ q ) Update (ṽ,σ p ,σ q ) ← (v, σ p , σ q ) end for Return (v, σ p , σ q ) and violation cost φ(σ p , σ q ) end for must discretize the homotopy path ρ(λ) = 0 and solve a finite number of homotopy-OPF problems. Let this discretization of homotopy path be represented by the set Λ : Then, the first homotopy-OPF problem, H λ o , is initialized as the solution to the base-OPF problem. The series of homotopy-OPF problems are solved via local search methods, where the solution to the previous homotopy-OPF problem is utilized as the initial point; see Algorithms 1 and 2 for complete details.
In this paper, we assume that the base-OPF has a unique global minimum that is available (known). The availability of a global minimum is a reasonable assumption because a good initial point is usually provided for the base-OPF, and also because more time is allocated to solving it compared to a large number of contingency-OPF problems for different outages, allowing the use of various convex relaxation techniques for the base-OPF. Now, note that the global minimum of the base-OPF is also a global minimum of H λ o because at λ = λ o , the parameters of the problem corresponds to the pre-contingency network.
If the optimal violation cost for H λ o is nonzero, the global minimum will be unique with overwhelming probability. Furthermore, even if the violation cost is zero, it will immediately become nonzero during the next homotopy iteration if removing a line or generator introduces inflexibilities that the network cannot accommodate. In fact, these near-infeasible problems where a contingency will make the system "stressed" are the cases where homotopy can be useful and are the focus of this paper. Later in the paper, we will present a rigorous result showing that the uniqueness of the global minimum is a generic property for H λ .
To develop intuition on when a homotopy method may or may not lead to the global solution, we consider the basin of attraction of the global solution to the contingency-OPF problem. The "basin of attraction" of a local solution is the set of initial points that lead to the solution using an iterative search method, and its size is dependent on both the problem geometry and the choice of algorithm.
At each iteration of homotopy, the problem is initialized at the solution to the previous problem with the hope that the previous point will be in the Under the uniqueness assumption mentioned above, the solution to some intermediate problem will enter and remain within the basin of attraction of the global solution to the final problem, and we will obtain the global minimum of the final problem. In the next section, we develop rigorous theoretical results for this idea.

IV. ANALYSIS OF HOMOTOPY PATHS
While probability-one homotopy methods almost surely guarantee algorithm convergence, they do not necessarily result in convergence to a global minimum [13]. In Section I-A, we offered two examples of nonconvex optimization: one in which the homotopy method resulted in the global minimum and another in which the homotopy method resulted in a non-global local minimum (see Figure 1). In this section, we describe a theoretical framework that describes when homotopy can be used to obtain a global minimum. We apply this framework to analyze the performance of homotopy-OPF in finding the global solution of the contingency-OPF. The results developed in this section have implications for homotopy methods in a broad range of optimization problems.
If the global solution of H λ is unique and satisfies strong secondorder conditions (SSOC) for every λ on the homotopy path ρ(λ) = 0, then a sufficiently small ∆λ will ensure that Algorithms 1 and 2 will find the global solution to H λ for all λ ∈ Λ.
The strong second-order conditions are similar to the second-order sufficient conditions for local optimality but with the addition of the strict complementary slackness condition and the linear independence of the active constraints [28]. The strict complementary slackness condition implies the lack of stationary points on the boundary of the feasible set. Furthermore, the SSOC is known to guarantee that a point is a strong local minimzer. Remark 1. We make the assumption that H λ along ρ(λ) has a unique global solution satisfying SSOC. In the next subsection, we argue that this is a reasonable assumption to make. In addition, this assumption on the global solution can be replaced by the "connectivity" of the set of all global solutions (this allows having infinitely many possible solutions for post-contingency OPF with zero violation cost).
In what follows, we show that the uniqueness of the global minimum is a generic property for H λ .

A. Genericity of Unique Global Minimizer with SSOC
Recall that a set S ⊂ R n has (Lebesgue) measure zero if for every > 0, S can be covered by a countable union of n-cubes, the sum of whose measures is less than . A property that holds except on a subset whose Lebesgue measure is zero is said to be satisfied generically or hold for almost all. In this subsection, we will show that the Homotopy-OPF generically has a unique global minimizer that satisfies SSOC .
Consider the following family of problems, which adds a linear perturbation to the objective of the Homotopy-OPF: where f : R 5|V| ×R → R, h : R 5|V| ×R → R 2|V| are smooth functions and the parameters (λ, ω) belong to an open set U ⊂ R × R 5|V| . We call this problem the extended homotopy-OPF. Here, represents the dimension of the parameter λ, which can be equal to either 4|E| (for line contingencies) or 2|V| (for generator contingencies). Then, using the results from [28], we can easily derive the following lemma: Lemma 1. Suppose the following two conditions are satisfied: 1) The function λ → h(x, λ) is of full rank 2|V| for all x ∈ Ψ at every λ 1 .
2) The set Ψ is a cyrtohedron and the set U is an open set. Then, for almost all (λ, ω) except those in a set U ⊂ U of measure zero, H λ,ω has a unique global minimizer satisfying SSOC. In fact, for every (λ, ω) ∈ U \ U , H λ,ω cannot achieve the same objective value at any two distinct critical points.
The concept of a cyrtohedra was first introduced in [29] and it captures a class of sets whose boundaries are a union of countably many smooth manifolds pieced together. A few main examples of cyrtohedra include polyhedral convex sets, submanifolds, submanifolds with boundaries, and manifolds with corners. In our case, the set Ψ is naturally a cyrtohedra and therefore we only have to verify the first condition. The next lemma proves that the condition can be easily verified for the line outage contingency.
where the column and row indices represent the lines and the nodes of the power system, respectively. If J has full column rank, then the function λ → h(x, λ) associated with the line outage homotopy method is of full rank 2|V|.
A similar result is proved for generator outage contingencies.

Lemma 3. Define the matrix
where both the column and row indices represent the nodes of the power system network. If M has full rank, then the function λ → h(x, λ) associated with the generator outage homotopy method is of full rank 2|V|.
The result implies that the first condition of Lemma 1 is satisfied if: (i) the pre-contingency real power generations and the post-contingency real power generations are different and (ii) the pre-contingency reactive power demands and the postcontingency reactive power demands are different. Note that this does not necessarily hold true because some real power generations are supposed to be fixed even after the contingency (same for reactive power demand). However, we can address this issue by allowing P g,f i (Q d,f i ) to take on a value within a small interval around P g,o i (Q d,o i ) whenever we want the two values to be close to each other.
Note that the linear perturbation term in H λ,ω is a mathematically necessary device that allows us to prove generic uniqueness of a family of nonlinear optimization problems. Ultimately, we will only consider very small perturbations so that H λ,ω closely resembles H λ . Hereby, letρ(λ, ω) = 0 represent an extended homotopy path that is continuous in the space of (λ, ω) that starts from the point (λ o , ω o ), which corresponds to the perturbed base-OPF, and ends at the point (λ f , ω f ), which corresponds to the perturbed contingency-OPF. Using the lemmas above, we arrive at the following main theoretical result of the paper:  Then, for every δ > 0, there exist sufficiently small > 0 and ∆λ > 0, and a global minimizer of H λ f , x * , such that solving the extended homotopy-OPF with Algorithm 1 (or 2) along the discretization ofρ(λ, ω) = 0 will find a pointx satisfying x − x * ≤ δ.
Here, solving the extended homotopy-OPF with Algorithm 1 (or 2) means sampling ω i randomly from the set B( ) at each instance of the discretized Λ i , and solving H Λ i ,ω i instead of H Λ i , which is the basic form presented in the algorithms. According to Theorem 4, the success of homotopy in finding the global optimum depends on the geometry of the set U and whether or not the homotopy path intersects with it. It is also shown that U is a measure zero set, which increases the likelihood of the existence of such non-intersecting path.

B. Geometry of the homotopy path: Two-bus example
In order to illustrate the previous ideas, we consider a simple homotopy-OPF example on a two-bus system. The line connecting the two buses has the admittance y = Gγ − jBβ, and there is a lower bound Q min on the reactive power injections at both buses. In this two-bus example, we consider the objective function: (σ p 1 ) 2 + c(σ p 2 ) 2 . Furthermore, assume that: where ∆ = tan −1 (Bβ/Gγ) and q(·) denotes the reactive power injection as a function of solely the angle difference, which is due to the fact that voltage magnitudes are fixed. Note that the second constraint on the angle difference is reasonable for the secure operation of power systems and is also used in [5] in order to restrict the two-bus real power injection region to be the Pareto front of the original feasible region. Geometrically, the feasible set of the two-bus injection region is the Pareto front of an ellipse, which is partially removed due to the reactive power constraints (the details can be found in [5]). Let P g,b i denote the real power generation at bus i obtained from the base-OPF solution. The following lemma characterizes the set of homotopy parameters for which there are at least two global solutions.
Lemma 5. Denote α = cos −1 −Q min +Bβ |y| , and define two polynomial functions of λ = (γ, β) as follows: Then, the set of parameters leading to multiple global minimizers, U , can be characterized as: · Ω 1 (γ, β) = 0} The set U in Lemma 5 for a particular instance of the example is depicted in Figure 2. As we can observe, U is a measure zero set in the two-dimensional parameter space, and it is possible to design an effective homotopy path. Note that the linear perturbation term in H λ,ω is a mathematical device used to prove generic uniqueness of the global minimizer for a family of problems. The characterization of U in Lemma 5 did not require the linear perturbation. However, this means that particular instances of the example may not lead to the result that we desire. For instance, if c = 1 and P g,b in the above example, U is no longer a measure-zero set.

C. Alternative characterization of desirable homotopy path
From the previous subsections, we learned that choosing a desirable homotopy path is directly correlated with the success of the proposed homotopy method. More specifically, the homotopy path must not intersect with a (measure zero) set consisting of parameters that lead to multiple global minima. Although the set being of measure zero helps, such nonintersecting path does not always exist. Therefore, in this subsection, we provide an additional perspective to analyze the desirable homotopy path using algebraic geometry.
Note that even though Algorithms 1 and 2 work on a discretized homotopy path, their analyses require working on the continuous path. In Figure 3, we present an example in which homotopy fails to find the global solution of the final problem. The major cause of this breakdown is the emergence of two global solutions in the (continuous) homotopy path, which is followed by a change in the relative positions of the global solution and next best local solution. In order to better characterize this, consider the KKT conditions for the homotopy-OPF problem defined in (3): where µ is the vector of dual variables. We assume that constraint qualification holds for the problem H(λ) defined in (3) for all λ ∈ S, which implies the KKT conditions are satisfied for every local minimum. The set S is defined in Theorem 4. For a given λ, let X (λ) be the set of all x that satisfy the KKT conditions in equation (6). Note that the goal is not to solve the KKT conditions directly but is merely to use them as a necessary condition for all local solutions. Before proceeding to the next theorem of this paper, we make one basic assumption on the KKT conditions and define a concept called the "dividing midpoint".
Assumption 1. The cardinality of set X (λ) as a function of λ is constant for all λ ∈ Λ. Also, SSOC holds for every point in X (λ) for all λ ∈ Λ.
The first part of Assumption 1 is essential to guarantee that a local solution is not suddenly created or disappeared along the homotopy path, in which case we cannot trace it back to the local solutions of the original problem to track it. Using techniques in algebraic geometry, one can study the satisfaction of the first part of Assumption 1 [30]. The second part of Assumption 1 is to make sure that we can apply the implicit function theorem. In [28], it is shown that this is generically true under mild regularity conditions.
Furthermore, let a be the midpoint objective value of the first and second best KKT points. In other words, Define U to be the set of all λ for which there exists a KKT point with the objective value equal to a: Here, we define a to be the "dividing midpoint" between f (x (1) , λ o ) and f (x (2) , λ o ). In practice, a wide range of values that are slightly above or below the point a would lead to the same implications. Theorem 6. Let ρ(λ) = 0 be a homotopy path of λ with two end-points λ o and λ f . If the set {λ | ρ(λ) = 0} does not intersect with the set U , then the homotopy problem H λ has a unique global minimum for all values of λ along the path ρ(λ) = 0.
Directly analyzing the geometry of the set U is cumbersome. Therefore, we introduce a method to certify whether a path is a successful homotopy path or not. The following theorem offers a dual certificate.
Note that the dual problem is convex and easy to obtain for certain problems, such as in the case where homotopy-OPF is cast as a non-convex quadratically-constrained quadratic program. In essence, finding a dual feasible point (ξ 1 , ξ 2 , ξ 3 , ξ 4 ) for which d(ξ 1 , ξ 2 , ξ 3 , ξ 4 ) > 0 provides a certificate that guarantees that the homotopy path ρ(λ) will never intersect with set U . Then, by Theorem 6, we can conclude that the homotopy method will have a unique global minimum along its path and therefore Algorithms 1 and 2 are able to solve contingency-OPF to global optimality using iterative local search due to Proposition 1. Note that the results of Theorem 6 are still valid if one adds valid inequalities to (P ) to increase the likelihood of the existence of a desirable dual feasible point via reducing the duality gap [31].

V. SIMULATIONS
Since theoretically proving that a given homotopy path works for a large-scale system is difficult, we numerically evaluate different homotopy paths and discretizations. In this section, we present simulations of different line and generator outage scenarios on different networks. From these simulations, we can evaluate when using homotopy to solve contingency-OPF is more or less useful than using a one-shot method.
In order to implement the contingency-OPF using the MAT-POWER format [32], we introduce virtual generators that model the violations of real and reactive power balances at all nodes after an outage occurs. These violations are penalized in a modified objective function. The benefit of this formulation is that there always exists a feasible solution to contingency-OPF. By adding power generation flexibility with virtual generators, we aim to find a feasible point (equivalent to a zero objective value) or an infeasible point for the In this case, all three homotopy schemes converge to a violation cost that is lower than the one obtained by the one-shot method. The bottom figure shows a 2 line outage (line out IDs: 17 and 342) in the 3120-bus Polish network (case3120sp) with real and reactive power demand scaled up by 10%. In this case, the one-shot and 3-and 5-iteration homotopy schemes converge to a violation cost that is lower than that obtained by the 10-iteration method.
network but with the minimum violations (such solutions could yet be implemented via corrective actions taken by realtime feedback controllers). To solve each of the homotopy simulations, we use the MATPOWER Interior Point Solver (MIPS) [33]. Figures 4 and 5 show the evolution of the violation cost over iterations of homotopy, compared to the violation cost of the one-shot method. Figure 4 shows line outage scenarios on the 3012-bus and 3120-bus Polish networks [32]. For these line outages, we implement a homotopy path that decreases γ ij , β ij simultaneously from 1 to 0 at each of the outed lines (i, j) ∈ E. We look at three possible discretizations of this homotopy path, using 3, 5, and 10 iterations of homotopy and compare them to the one-shot method. Figure 5 shows generator outage scenarios on the 89-bus and 1354-bus PEGASE networks [34], [35]. For these generator outages, we implement a homotopy path which decreases λ from [1 |V| , 1 |V| ] to [0 |V| , 0 |V| ] uniformly throughout the iterations. For this homotopy path, we also consider 3, 5, and 10 iterations of homotopy. From these figures, we can see that the final violation cost obtained using the given homotopy paths can vary significantly depending on the number of iterations (i.e. ∆λ) of homotopy-OPF. See the Appendix for more simulations, including simulations that compare different homotopy paths rather than discretizations.
In order to formally compare the performance of homotopy versus the one-shot method, we say that homotopy "outperforms" the one-shot method if either of the following are true: 1) If the homotopy scheme converges and the one-shot method does not converge. 2) If the homotopy scheme converges to a value that is better than that of the one-shot method by at least 0.01% of the optimal base-OPF cost. For the 1354-bus PEGASE network, we tested 1, 2, and 3 line and generator outages, testing 100 simulations of each type of outage. The homotopy paths for these line and  generator outages are the same as those described for the simulations in Figures 4 and 5. The percent of simulations where homotopy outperformed the one-shot method is given in Table I for the network with base-level demand, and in Table II for the network with demand scaled up by 10%. It can be observed that for the line outage contingencies, the homotopy methods appear to be more useful when the demand is higher. This is likely because the increased demand makes the problem harder, and thus homotopy is more useful. However, the inverse appears true for the generator outage scenarios, i.e. the homotopy methods appear to be more useful when demand is at the base-level. This could be because the removal of a generator could lead to many possibilities for operating the post-contingency network in a lower demand scenario, which may introduce bad local minima.

VI. CONCLUSIONS
This paper studies the contingency-OPF problem, which is used to find an optimal operating point in the case of a line or generator outage. Unlike the base-OPF problem that is a single optimization problem, there are many contingency-OPF problems that should all be solved in a short period of time. Recognizing that the contingency-OPF problem is a more challenging version of the classical OPF problem, we introduce a new homotopy method to find the best solution of the contingency-OPF problem. This method involves solving a series of intermediate homotopy-OPF problems using simple local search methods, and we study conditions that guarantee convergence to a global solution of the contingency-OPF. We perform simulations on real-world networks and show that the proposed homotopy method can result in a lower value of the objective.

A. Proof of Proposition 1
Let x * 1 denote the unique global solution satisfying SSOC for the problem H Λ 1 (also, recall that Λ 1 = λ o ). Using an argument relying on the implicit function theorem [36], it follows that for each (x * 1 , Λ 1 ) pair, there exists a neighborhood U 1 around Λ 1 and a neighborhood X 1 around x * 1 , and there is a differentiable function x 1 (λ) defined for λ ∈ U 1 such that 1) 2) For each λ ∈ U 1 , x 1 (λ) is the unique point in X 1 satisfying the SSOC for H λ . Now, suppose that ∆λ is small enough so that Λ 2 ∈ U 1 . Then, since x 1 (λ) is a continuous function and there is no λ on the path ρ(λ) = 0 such that H λ has two global minimizers, x 1 (Λ 2 ) becomes the unique global minimizer satisfying SSOC for the next OPF problem, H Λ 2 . The same logic can be applied for all Λ i , and by induction we have proved the result.

B. Proof of Lemma 1
By Proposition 4 of [28], the family of optimization problems has a unique global minimizer satisfying SSOC for all parameters (λ, ω) ∈ U except on a set of measure zero if (i) for all x 1 = x 2 , and for all ω, the function ω → f (x 1 , λ, ω)−f (x 2 , λ, ω) is of rank one at all λ, (ii) the function λ → h(x, λ) is of full rank 2|V| for all x at every ω, and (iii) the fixed set Ψ is a cyrtohedron and U is an open set. It is straightforward to check that iff (x, λ, ω) = f (x, λ) + ω T x, condition (i) is satisfied. Conditions (ii) and (iii) are given as assumptions, which completes the proof.

C. Proof of Lemma 2
The rank of the function λ → h(x, λ) is the rank of its Jacobian (w.r.t. λ). Therefore, we analyze the Jacobian of h(x, λ) = h(x, [γ, β, γ sh , β sh ]) with respect to [γ, β, γ sh , β sh ]. From Section II-B and III-A, we know that h consists of two types of functions, h 1 and h 2 (corresponding to the real power flow equations and the reactive power flow equations, respectively), whose i-th elements are defined by: We focus on the submatrix of the Jacobian that consists only of the derivatives of h 1 i and h 2 i with respect to γ sh and β sh (denote this as J). This is because if this submatrix has full column rank, then the full Jacobian also has full column rank. First, we notice that the Jacobian of h 1 with respect to β sh and the Jacobian of h 2 with respect to γ sh are equal to zero. Therefore, J can be expressed as a 2 × 2 block matrix of the form J = J 1 0 0 J 2 ∈ R 2|E|×2|V| , where J 1 corresponds to the Jacobian of h 1 with respect to γ sh and J 2 corresponds to the Jacobian of h 2 with respect to β sh .
For line outage contingencies, γ sh and β sh are parameters indexed by the line number. Hereby, let J 1 ((i,j),k) refer to the element of J 1 that is located at the (i, j)-th row and the k-th column (the row index representing the line and the column index representing the bus number). For example, J 1 denotes the partial derivative of the real power flow equation at bus k with respect to the shunt susceptance parameter at line (i, j). The same goes for J 2 .
Then, directly from basic calculus, we can derive the following form for the matrix J: Therefore, if J has full column rank, so will the Jacobian of the function λ → h(x, λ), which completes the proof.

D. Proof of Lemma 3
Similar to the proof of Lemma 2, we analyze the Jacobian of h(x, λ) = h(x, [γ, β]) with respect to [γ, β]. From Section II-B and III-A, we know that h consists of two types of functions, h 1 and h 2 (corresponding to the real power flow equations and the reactive power flow equations, respectively), whose i-th elements are defined by: First, we notice that the Jacobian of h 1 with respect to β and the Jacobian of h 2 with respect to γ are equal to zero. Therefore, M can be expressed as a 2 × 2 block matrix corresponds to the Jacobian of h 1 with respect to γ and M 2 corresponds to the Jacobian of h 2 with respect to β. For generator outage contingencies, γ and β are parameters indexed by the bus number (because we assume each bus has exactly one generator). Hereby, let M 1 i,j refer to the element of M 1 that is located at the i-th row and the j-th column. In other words, M 1 i,j denotes the partial derivative of the real power flow equation at bus i with respect to the γ parameter at bus j. The same goes for M 2 .
Then, directly from basic calculus, we can derive the following form for the matrix M : Therefore, if M has full column rank, so will the Jacobian of the function λ → h(x, λ), which completes the proof.

E. Proof of Theorem 4
The first statement on H λ,ω having a unique global minimizer satisfying SSOC follows directly from applying Lemma 1. The functions λ → h(x, λ) is of full rank 2|V| due to Lemmas 2 and 3, and this in turn satisfies the first condition of Lemma 1. As discussed in Section IV, the set Ψ is a cyrtohedron, and the set U is defined to be an open set for any > 0 by the assumptions of this theorem. In other words, the second condition of Lemma 1 is also satisfied. Therefore, we can conclude that for any value of > 0, H λ,ω has a unique global minimizer satisfying SSOC for every (λ, ω) ∈ U \ U where U ⊂ U is of measure zero.
Next, we move on to prove that solving H λ,ω where ω belongs to a sufficiently small open ball of radius will be able to track the global minimizers of H λ . In order to do so, we first mention an existing result which states that small perturbations on the function only change the local minima by a small amount. Lemma 8. [14] Let f : R n → R and g : R n → R be twice continuously differentiable functions, and let x * be a non-singular local minimum of f. Then, there exist > 0 and δ > 0 such that for all¯ ∈ [0, ) the function f (x) +¯ g(x) has a unique local minimum x * * within the open ball {x | x − x * < δ} Let x * andx * denote global minima of H λ and H λ,ω , respectively. We note that the problem H λ,ω simply represents a linear perturbation on the objective function of H λ . The feasible sets of the two problems are exactly the same. The objective function of H λ,ω can be expressed as: where ω denotes the 2 norm of ω. Since the set {(λ, ω) | ρ(λ, ω) = 0} does not intersect with U (except for possibly at λ = λ f ),x * is the unique global minimum of H λ,ω satisfying SSOC. A local minimum that satisfies SSOC is non-singular. Therefore, by Lemma 8, we know that there exist and δ such that for all ω ≤ , the function f (x, λ) has a unique local minimum within the open ball B(x * , δ) = {x | x −x * < δ}. Note that we can use Lemma 8 (a result on an unconstrained optimization problem) because SSOC conditions (specifically, strict complementary slackness) hold atx * . In fact, we can make small enough so that the unique local minimum of f (x, λ) within B(x * , δ) is also a global minimum and satisfies SSOC. Finally, similar to the proof of proposition 1, we can show that the sequence of problems starting from H Λ 1 ,ω 1 can track the unique global minimizers of the intermediate problems along ρ(λ, ω) = 0 all the way through H Λ T ,ω T , given a small enough ∆λ along with the small enough . Combining this with the statements made above, we conclude that there exist sufficiently small and ∆λ such that solving H λ,ω along the discretization of ρ(λ, ω) = 0 will find a point that is δ − close to the global minimizer of H λ .

F. Proof of Lemma 5
Let us start with the equation for the reactive power injections. Let θ 1 and θ 2 denote the voltage phasor angles at buses 1 and 2, respectively. Let the real and reactive power injections at bus i be denoted by p inj i and q inj i , respectively. In this two-bus example, we consider the objective function: (σ p 1 ) 2 + c(σ p 2 ) 2 . Then after denoting θ = θ 1 − θ 2 , we have the following: A lower bound of Q min on q inj 1 results in the following: Q min ≤ Bβ − Gγ · sin θ − Bβ · cos θ Then, after rearranging and using trigonometry, we arrive at −Q min + Bβ ≥ Gγ · sin θ + Bβ · cos θ = (Gγ) 2 + (Bβ) 2 · cos (θ − ∆) where ∆ = tan −1 Gγ Bβ . and (p inj 1 ,p inj 2 ) are both globally optimal, their objective values must be equal. In other words, Rearranging the terms leads to 2 )(p inj 1 −p inj 1 ) = 0 Finally, substituting the definition of Ω 1 and Ω 2 , we arrive at (1 − c) · Ω 1 (γ, β) · Ω 2 (γ, β) − 2(P g, b 1 − P d 1 ) · Ω 1 (γ, β) + 2c(P g, b 2 − P d 2 ) · Ω 1 (γ, β) = 0 , we see a case where the 2, 5, and 10iteration homotopy methods result in an objective value much lower than that obtained by the one-shot method. For this scenario, by introducing even a 2-iteration homotopy scheme we outperform the one-shot method.   9. Performance of proposed homotopy method on the 3120-bus Polish network (case3120sp with original real and reactive power demand) with a single line outage (line out ID: 1602). Homotopy schemes 1 through 3 are tested with 10 iterations. In this example, homotopy schemes 1 and 2 result in an objective value higher than that obtained by the one-shot method, while the third homotopy scheme outperforms the one-shot method.