{'fraction_non_alphanumeric': 0.03574816220061655, 'fraction_numerical': 0.008032961821199905, 'mean_word_length': 4.827776818103299, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Incremental Model Building Homotopy Approach for Solving Exact AC-Constrained Optimal Power Flow

Alternating-Current Optimal Power Flow (AC-OPF) is framed as a NP-hard non-convex optimization problem that solves for the most economical dispatch of grid generation given the AC-network and device constraints. Although there are no standard methodologies for obtaining the global optimum for the problem, there is considerable interest from planning and operational engineers in finding a local optimum. Nonetheless, solving for the local optima of a large AC-OPF problem is challenging and time-intensive, as none of the leading non-linear optimization toolboxes can provide any timely guarantees of convergence. To provide robust local convergence for large complex systems, we introduce a homotopy-based approach that solves a sequence of primal-dual interior point problems. We utilize the physics of the grid to develop the proposed homotopy method and demonstrate the efficacy of this approach on U.S. Eastern Interconnection sized test networks.


Introduction
Alternating Current Optimal Power Flow (AC-OPF) attempts to solve for the most economical dispatch of the grid while satisfying AC network constraints and enforcing device limits. Posed as a non-convex, nonlinear optimization problem , AC-OPF is NP-hard with no methodology for guaranteeing a global optimum in polynomial time for general networks even though some work [1]- [2] have found the global optimum for a restricted set of small networks.
Instead, many methodologies have been proposed and developed to solve for the local optimum or an approximation of the AC-OPF problem due to increasing demand from grid planners and operators [3]. These methodologies can be divided into three main classes: i) those that solve for the local optimum of the original AC-OPF problem [4]- [5]; ii) those that solve the linear approximation of the AC-OPF problem (e.g. Decoupled-OPF) [6]- [8]; and iii) those that solve the convex relaxation or restriction of the AC-OPF problem [9]- [12].
Of these three classes of methodologies, in this paper, we focus on obtaining the local optimum solutions for the AC-OPF problems primarily because the local optimum ensures a feasible dispatch with AC network constraints satisfied and device limits enforced. In planning and operation, security and reliability are considered more critical than the cost of the electricity and obtaining a sub-optimal feasible solution to AC-OPF is generally preferred to an AC-infeasible solution (due to relaxations) that is more economical.
Due to these factors, there has been a recent surge in methods that locate local optimum solutions for the AC-OPF problem [4]- [5]. In particular, a large-scale effort was driven by an Advanced Research Project Agency -Energy (ARPA-E) grid optimization challenge [13]. The goal of the competition was to locate the most economical dispatch of the AC-OPF while ensuring feasibility across a set of contingencies (also known as Security Constrained Optimal Power Flow problem). Many state-of-the-art methods for obtaining local optimum solutions for non-convex optimization problems were presented in the competition [13], with the most successful ones based on standard nonlinear optimization toolboxes such as the Matpower Interior Point Solver (MIPS) [14] and Interior Point Optimizer (IPOPT) [15]. The approaches utilized power-system specific knowledge to achieve robust convergence. An alternative approach made use of a circuit-based formulation with currents and voltages as state variables [5] that applied circuit-simulation methods to achieve robust convergence. Almost all of these approaches demonstrated the ability to solve real-size AC-OPF problems within the SC-OPF problem; nonetheless, the results of the challenge [13] suggest that there is value in further exploration. Most importantly, since the optimization framework also serves as a basis for planning engineers, a robust approach cannot always rely on starting with good initial conditions.
Various previous works have attempted to develop a robust AC-OPF approach by utilizing a homotopy method technique [16] that traces a path in the solution space by incrementally modifying the non-linear constraints from a trivial problem to the original AC-OPF. This method requires a continuous path to the solution with a valid solution at each increment. However, previous works [18]- [24] are unable to ensure a valid feasible solution at each increment and few create discontinuities in the homotopy path due to the use of DC-OPF as an initial solution [19]. As a result, many homotopy methods used for the power grid today [18]- [24] are unable to scale to large networks.
One homotopy method that is effective for largescale power-flow problems is Tx-stepping [25]. Txstepping is based on an approach that is analogous to Gmin stepping for circuits [27] and aims to short the series elements (e.g. lines, transformer, etc.) in the system by adding large parallel admittances at first to obtain a trivial solution. Throughout the homotopy path then, the added admittances are gradually removed to solve the original network. Here for grid optimization we employ a more comprehensive incremental formation of the final system starting from a trivial one that not only scales the transmission line and transformer parameters, but also shunts, generators, and loads and corresponding bounds. This effectively corresponds to an incremental formulation of the entire grid from a trivial one that mimics turning on the grid. Importantly, this incremental stepping must be performed such that there is a valid grid model at each step, which we ensure by using a slack injection-based technique that maintains feasibility throughout the homotopy method. We refer to this homotopy as IMB, Incremental Model Building.
We introduce a framework that implements the IMB approach within the Simulation of Unified Grid Analyses and Renewables (SUGAR) engine [5], [17]. We overcome convergence challenges that generalpurpose optimization toolboxes struggle with by utilizing the known properties of the physics-based grid models to ensure convergence, even when confronted with ill-suited initial conditions. Importantly, the proposed homotopy approach is generic and can be applied to achieve robust convergence for other grid optimization problems as well.
Section 3 covers the AC-OPF formulation along with primal-dual interior point (PDIP) solution methodology. Section 4 introduces the novel homotopy method and Section 5 describes the algorithmic methodology following the proposed approach. We conclude the paper with AC-OPF results for large and complex networks and compares it against state-of-theart standard optimization toolboxes in Section 6.

Nomenclature
, , coefficients of quadratic cost function for generator .
, , , Real and imaginary currents injections by a generator.
, , , Real and imaginary currents injections by a load.
, Real and imaginary voltage at node .
, Real and imaginary voltage across nodes and .
, Components of admittance matrix.
, Real and reactive power of generation .
, Real and reactive power of constant power demand .
Square of the voltage magnitude at bus .
Current flow in a series element .
Apparent power flow in a series element .
, Bounds on generator real power.
, Bounds on generator reactive power.
, Bounds on voltage magnitude at bus . Homotopy factor.
Scaling factor of admittance during homotopy.

Current Voltage Formulation
Traditionally, the non-linear AC-network constraints of an AC-OPF problem are represented by power mismatch at each bus; however, the trigonometric terms that model network line flows and network constraints introduces steep non-linearities in the solution space when scaling to large systems. A growing number of approaches [3], [5] are adopting a current-voltage formulation in which network constraints are expressed as Kirchhoff's Current Laws (KCL). The currentvoltage (I-V) models have better convergence properties for the AC-OPF problem than the traditional power mismatch-based formulations due to the reduced nonlinearities in the network constraints [3]. Unlike the traditional formulation where the network constraints are non-linear with trigonometric terms, the I-V formulation has linear network constraints with nonlinearities stemming from the injection models (loads and generation). Therefore, to minimize non-linearities in the formulation, we will utilize the I-V representation of the AC-OPF problem.
The optimization formulation for the AC-OPF problem of a power grid network Ɲ consisting of a set of generators and load demands connected to a set of buses or nodes in the grid is given by (1). The nodes in the system are connected by a set of network elements, and . The objective function for AC-OPF of minimizing generation cost is given by ℱ ( ), defined by sum of quadratic functions given by a set of coefficients { , , }. The equality constraints represent the network constraints given by (1b)-(1h). Inequality constraints represent physical bounds on the devices given by (1i)-(1l). Some bounds include system limits that are based on grid stability; for example flow constraints can be both thermal or stability based (see 1l) and are further discussed in Section 4.4. In the problem formulation, system topology parameters, cost, and device and voltage bounds are given parameters whereas generator real and reactive power output along with voltage-setpoint are the decision variables ( ⊆ = { , , }) over whose range the problem is optimized. More detailed formulations can also include transformer tap and phase-shifters ( , ) along with shunt positions as decision variables.
subject to: With the non-linear optimization formulation devised, we can construct a framework for solving for a local minimum using Primal-Dual Interior Point method.

Solution Framework for AC-OPF Problem
To solve the non-convex optimization problem given in (1), we implement the primal-dual interior point (PDIP) approach (see [26] for details). The PDIP approach is the basis for further heuristics developed within this paper that enables robust and scalable convergence. PDIP algorithms apply the search direction from Newton's method to iteratively solve a set perturbed (Karush Kuhn Tucker) KKT conditions. The perturbed KKT conditions represent a relaxation of first-order optimality conditions that are necessary to obtain a local optimum for the AC-OPF problem. We begin by formulating the Lagrangian: where, ( ) represents the vector of equality constraints given by (1b)-(1h) (in case of AC-OPF) and ℎ( ) represents the vector of inequality constraints (1i)-(1l). represents the vector of primary variables whereas and represent the vector of dual and slack variables corresponding to equality and inequality constraints, respectively. In the case of AC-OPF, is a vector of grid states that include real and reactive power output of the generators ( , ), vector of real and imaginary components of grid voltages ( , ), vector of line and transformer current and power flows ( , ) and a vector of generator set-points ( ). The set of perturbed KKT conditions (perturbed due to the relaxation in complementarity slackness constraint) corresponding to (2) can now be written as: where ∘ represents element-wise multiplication of the vector elements and represents a vector of complementarity slackness tolerance.
To obtain a stationary point ( * , * , * ) that satisfies (3), PDIP linearizes and iteratively solves the equations corresponding to (3a)-(3d) using Newton's method. Between two NR iterates, primal and dual feasibility given by (3e)-(3f) is satisfied through heuristics. The linearized matrix at iteration can be written as follows: where ≔ diag( ), ≔ diag( ) and = ∇ 2 ℒ( , , ). In general, the linearized matrix is reduced and 's eliminated to solve lower dimension matrix within the inner loop of the Newton's method. A descent direction is guaranteed by ensuring that the top left block (given by ) of the iteration matrix is positive definite [15] and step-length of the iterate is often controlled via backtracking line search [26].
While this straightforward PDIP solution methodology could be used to solve small networks, it is unlikely to be effective for complex large networks with steep nonlinearities. Therefore, to ensure convergence for such realistic networks, heuristics are needed. Most commercial tools have such heuristics embedded within them. For instance, MATLAB optimization toolbox makes use of interior trust region methods [32], whereas IPOPT makes use of filter line search-based methods along with second order correction terms [15]. Nonetheless, these methods also diverge for some of the hardest networks as these do not intrinsically utilize grid physics within their solvers. Therefore, to ensure robust convergence of hardest AC-OPF networks, we augment the naïve PDIP algorithm with a novel IMB homotopy method that utilizes the physics of the electric grid to robustly obtain a local optimum.

Homotopy Methods
Homotopy methods [16] are commonly used to solve complex large-scale non-convex optimization problems pertaining to many applications. However, so far, the use of such methods for solving the power grid optimization problems is limited [17]- [24]. These methods embed a scalar homotopy factor, , into the problem in an effort to relax the non-linearities. Varying the homotopy factor, : 1 → 0, replaces the original problem with a set of sub-problems represented by each increment of the homotopy factor and are sequentially solved with the following properties: i) the first sub-problem has a trivial solution and ii) each subsequent problem has a solution very close to the solution of the prior sub-problem wherein the final subproblem represents the original problem. This second property exploits the Newton-Raphson quadratic convergence properties thereby allowing faster convergence. Mathematically this can be described via the following expression: where ∈ [0, 1]. The method begins by embedding the homotopy factor into the original problem, effectively replacing the original problem ℱ( ) = 0 with ℋ( , ) = 0. The equation set ( ) is a representation of the system that has a trivial solution. The homotopy factor has the value of 1 for the first sub-problem that corresponds to a trivial problem ( ). Iteratively the homotopy factor, , is reduced to 0 which then represents the original problem ℱ( ). In the homotopy path between the trivial problem ( ) and the original problem ℱ( ) lies a sequence of sub-problems that traces a path in the solution space. Importantly, the homotopy method requires a valid solution at each step of the homotopy path to ensure convergence.
Previously, the homotopy method has been proven to be an effective tool to achieve robust convergence for complex power flow analysis problems in the form of a Tx-stepping in [25]. However, power flow analysis assumes a fixed generation and solely applying Txstepping in an AC-OPF setting proves ineffective as generation is variable in the problem definition. In this paper, we extend the Tx-stepping homotopy method to develop a broader more comprehensive homotopy method to solve optimization problems. Specifically, we target robust solution of AC-OPF problems using our Incremental Model Building (IMB) approach.

Tx-Stepping Homotopy Method
In previous work [25] we developed a homotopy method, Tx stepping, that solved for any positive sequence or three-phase power flow (PF) problem independent of the given initial conditions. Tx-stepping is based on embedding a homotopy factor in the linear transmission network to virtually short the grid. This homotopy approach begins by solving an almost shorted system and then it gradually decreases the homotopy factor to obtain the solution to the original problem. We briefly describe how grid models are modified based on Tx-stepping method.

Transmission line models.
In Tx-stepping, the series elements in the system (transmission lines, transformers etc.) are "virtually" shorted by adding a large conductance ( ) and a large susceptance ( ) in parallel to each transmission line and transformer model in the system parameterized by the homotopy factor: where, is the set of all transformers and is the set of all the transmission lines in the system. and are the original line and transformer impedances and the � and � are the parameterized admittances during Tx-Stepping that are used while iterating from trivial problem to the original problem. The parameter is used as a scaling factor for the conductances ( ) and susceptances ( ). Evidently, the first homotopy factor will increase the conductance and susceptance of all the branches and transformers �( � � ), ∀ ∈ { , }� thereby almost shorting the system and incrementally the homotopy factor will decrease these values to the original network state at = 0.

Transformer phase shifters and taps.
A "virtually short" power system (at a homotopy factor of 1) will drive all the voltage magnitude and angle of all the buses to the same value driven by the reference bus. To ensure transformer taps and phase shifter angles are compliant, their turns ratios, ̂ and phase shift angles, � correspond to a magnitude of 1 pu and 0°, respectively. Subsequently, the homotopy factor is reduced to constrain the transformer tap and phase shifters to their original settings. This can be mathematically expressed as follows: Note that with following changes to the transmission line and transformer models, we have shown to solve any power flow problem robustly, independent of the choice of the initial conditions [17]. However, this approach is insufficient and unlikely to work for AC-OPF as the generation is variable and the physical limits of devices must be enforced. Therefore, to overcome these challenges, we devise a new homotopy method as an augmentation of Tx-stepping method.

Homotopy Models for Optimal Power Flow
Incremental Model Building (IMB) builds upon Txstepping to provide a homotopy method that is designed for optimizing power grids. The goal of the IMB homotopy method is to mimic gradually turning on the grid from a shorted no-load network. Specifically, we first scale down all the loads in the network to virtually zero in the first homotopy step by embedding a homotopy factor in load and other AC-OPF based grid models, thereby creating a viable solution model that satisfies the physics constraints. From there the embedded homotopy factor iteratively traces a path from the initial sub-problem, depicted in Figure 1 (left) to the final one as shown later in Figure 1 (right). Along with virtually shorting the transmission network as is done in Tx-stepping, the IMB approach embeds the homotopy factor into the load, generator, and shunt models to trace the solution path.

Load Models.
To mimic turning off the grid at first, we embed a homotopy factor ( ) that reduces a load's active and reactive powers to virtually zero at first and then gradually scales it back up as it traverses through the homotopy path, as shown by (10)- (11). The initial trivial problem with ̅ = min(0.0001, ), modifies the network to have virtually no load.

Generator Limits.
During the initial step with the "turned-off grid" homotopy step, the output of generator's real power is very close to zero ( ≅ 0, ∀ ∈ ). To ensure that this output is feasible and within the generator real power bounds, we adjust the upper and lower limits of the generator real and reactive power based on the homotopy parameter. This ensures feasible generator operation in the homotopy path even when the original lower bound of the generator real power is much greater than 0 ( ≫ 0) or the original higher bound is much lower than 0 (0 ≫ ).
where is a fixed constant (in our experience a value of 0.5 works well).

Shunt models.
Following the IMB technique, shunt devices are initially turned-off and then gradually turned on while traversing through the homotopy path to have its full capacity when solving the original problem: ( � ℎ , � ℎ ) are homotopy parameterized shunt impedances that replace the original impedances ( ℎ , ℎ ) at bus . Iteratively with the homotopy factor, the shunt admittances are restored to their original value.

Circuit-based Flow Constraints
Power and current-based limits are two common ways in which flow constraints for transmission line and transformers are represented within AC-OPF problem. They either represent thermal or stability limits for the series element. Commonly, current bounds are used for lines and power bounds are used for transformers. We first describe how we implement current bounds. Power bounds can be implemented similarly. Modeling the line and transformer current bounds require computing the real and imaginary currents through any given transmission line or transformer. To calculate the currents as in the case of circuit simulation, we append a voltage source with a voltage of 0 pu (representing an ammeter) in series with the series element, as shown by Figure 2. This results in addition of two unknown states to the set of equations, which corresponds to measured real and imaginary currents ( and ). Importantly, the voltage measurement source is a linear element that is able to measure the real and imaginary currents, of and without introducing any non-linearities. Now that these currents are state variables, we can add an additional equation (17) to calculate the square of magnitude of the current through the series element, ( ) 2 , where ∀ ∈ { , }.
Inserting (17) into the Lagrange function (2) with an associated dual variable, , , we can bound the current magnitude to a value ���� by adding following constraint: Similar methodology can be followed if transformer bounds are given in apparent power units, with an added equation for power calculated from receiving node and sending node: and a corresponding set of dual variables ( , , , ) and dual equation as (18).  Existing homotopy methods for power grids fail to ensure this criterion [19], [21]- [22]. To satisfy this criterion in our proposed homotopy method, we draw inspiration from the physics of the grid. In the AC-OPF problem formulation (see Section 3), the equality constraints represent the conservation of charge at each node in the grid (modeled through KCL) and the inequality constraints enforce the device limits.

Feasible Space
Infeasible Space

Without Slack Injections
Generally, it is possible that there may exist a PDIP subproblem corresponding to a homotopy factor in homotopy path ( ( )), whose solution may be "infeasible" when the required current cannot be supplied by the system to satisfy KCL while ensuring that the devices are within their physical limits. If such a scenario occurs during the homotopy path ( ), the homotopy method would fail as no feasible solution exists for a given sub-problem k (see Figure 3 bottom).
To overcome this challenge, we design the IMB homotopy method to guarantee a feasible solution throughout the homotopy path by adding slack injections ( , ) (also see [28]) to each node in the network that are multiplied by the homotopy factor ( , and , ). These represent variable current sources at each node in the network to satisfy any violation of KCL. The magnitudes are minimized in the objective function with a large weight to sufficiently penalize any slack currents. The KCL equations at each node in the network given by (1b)-(1c) in the AC-OPF formulation are modified during homotopy to include the slack currents as follows: By adding this homotopy factor based slack injections, we can vary the magnitudes of slack injections to always have a feasible solution for each step in the homotopy path. Equivalently, a feasible network implies that the Jacobian ( ℋ ′ ( )) is full rank at each point on homotopy curve ( ) as KCL is satisfied at each node due to the existence of slack injection sources. We further ensure that these injections are removed in the final sub-problem ( = 0) to solve the original network, as they are parameterized by the homotopy factor . Note that this assumes there exists at least one feasible solution to the original problem ( * , * , * ) that satisfies (3a) through (3d). In case there is a problem that does not have a feasible solution (one that satisfies (1b)-(1l) at = 0), then one can obtain its physically infeasible solution by satisfying only relaxed form of (1b)-(1c) by replacing , , in (20)- (21) with , where is a very small number.

Limiting for Primal-Dual Feasibility
For each homotopy step increment, we solve a PDIP sub-problem that must satisfy (3a) through (3f). In the IMB approach, we satisfy (3a)-(3d) through solving the set of non-linear equations following the damped Newton's method. We then satisfy (3e)-(3f) during each NR step through the technique described below. Primal feasibility (3e) is ensured when the primal variables, , do not violate their appropriate device limits, ℎ( ) ≤ 0. Similarly, the dual feasibility (3f) is maintained when the dual variables are non-negative, ≥ 0. To achieve the following, we employ a limiting heuristic similar to that of diodes in circuit simulation [27]. The limiting method ensures that any limited primal variable does not violate its limits (satisfies (3e)) during each NR iteration. To do so, every variable in the vector , has an associated damping factor, , ∈ (0,1] that modifies the updated variable, +1 at the th iteration. where +1 is the value at the next iteration that is determined by the step, Δ . To properly limit the variables during each iteration to satisfy the inequality constraint ( < ̅ ), the damping factor is calculated as where is a constant (typically of 0.95-0.99) to ensure that +1 does not hit its limit exactly. Importantly, each variable being limited calculates its own separate damping factor. Similarly, corresponding dual variables are limited to remain non-negative within each NR step: and vector of damping factor ( ) is calculated following this logic:

Convergence Notes on Homotopy Methods
Global convergence of homotopy methods require following conditions to be met [16]: i. Defined path for the homotopy method i.e. ( ) ∈ ℋ −1 (0) with ( , ) ∈ range( ) must be smooth and should exist. ii.
If a curve exists, then it should intersect the final solution at = 0. In the proposed IMB homotopy method, the first condition can be met through implicit function theorem and requires that the Jacobian (ℋ ′ ( )) of the homotopy function is of full rank for all values of along the curve. For AC-OPF, we ensure the full-rank of the iteration matrix by following the methodology in Section 4.5.
The second condition is more easily met and is linked to existence theorems in non-linear analyses [16]. If some boundary condition exists that prevents the curve from extending to infinity prior to intersecting the solution at = 0, then this condition is met. In our formulation, different limiting methods ensure that the solution at any point on the curve does not diverge and extend to infinity given a feasible solution exist for the defined problem.

SUGAR Optimization Algorithm with IMB Homotopy
Using the IMB homotopy algorithm with convergence guarantees under the hold of certain conditions, we integrated the method within the circuitbased solver power flow solver, SUGAR [17] to solve the AC-OPF problem. The framework begins by accepting an input file(s) to parse the network data, generator costs and device operating limits. Following this step, all system states are initialized based on the values in the input file or based on other heuristics and any out-of-bound states are reset to feasible values within the bounds. Next, the device models are generated for elements in the network file and a homotopy parameter is embedded within the models (see, [17] and [5] for details). For non-linear models, the non-linear Hessian terms are initialized based on the initial states. The perturbed KKT conditions corresponding to the grid optimization problem at initial homotopy factor (υ=1) are iteratively solved with Newton method, thereby providing the search direction to convergence. The homotopy factor is updated and the process is repeated with the final solution from the previous homotopy factor step used as the initial condition for the new sub-problem with updated homotopy factor. Within each NR iteration, the residual and error corresponding to the primal and dual constraints given in (3a)-(3d) are used as input heuristics to determine the homotopy factor step and the Newton step size. Also limiting heuristic based on Section 4.6 is used to ensure primal and dual feasibility given by (3e)-(3f). The algorithm is terminated once convergence is achieved for the primal and dual constraints at homotopy factor of 0 and primal and dual variables are in feasible space. ALGORITHM 1: procedure: 2: parse network data, generator cost, and device limits 3: initialize system states to be in the feasible region; assign tolerance 4: modify the device models and bounds to initialize homotopy (υ=1) 5: create the matrix structure 6: loop through all the linear devices and create linear admittance matrix 7: loop through all the non-linear devices and initialize Hessian matrix of non-linear terms based on initial state while homotopy factor (υ≠0) 8: while not converged do: Solve linearized KKT with homotopy factor, υ, using NR 7: enforce primal and dual feasibility 8: check error and residual profile 9: if error or residual profile diverging: 10: adjust homotopy factor (υ) (decrease homotopy step) 11: end if 12: if change in υ: update homotopy elements 13: update non-linear devices 14: end while; update homotopy factor (increase the step) 15: end while Cases shaded in grey represent congested networks also superscripted using C 3. Superscripts 1, 2, and 3 represent three different objective functions. All the input .m files are uploaded to public git repo.

Results
To demonstrate the efficacy of our approach, we optimize the AC dispatch (with AC-OPF) for multiple test networks ranging from 9k+ nodes pegase network [30] to 80k+ node synthetic test networks for the U.S. interconnection [31]. We run two sets of analyses on each network with one corresponding to a set of noncongested scenarios (unshaded in the table) where none of the transmission line limits are binding and another corresponding to a set of congested scenarios (shaded in grey in the table) where some of the transmission line limits are binding. We run these two sets of scenarios (in total 15 cases) and compare the obtained results against those obtained from three other solvers with the Hessian supplied through a user-defined function: • fmincon [30]  The following changes to the default settings are made for the commercial solvers: i) Increase the maximum iteration count for the IPOPT solver from 250 to 750; ii) Represent maximum flow limits in terms of currents (I) instead of default setting of real power (P) for congested cases; and iii) For non-congested scenarios, use default settings for flow limits of real power (this is because convergence was much worse when flow limits in terms of currents were used).
All of the .m testcase files with generation cost data, line flow limits and other case information are made publicly available in a Gitlab repository [29]. Also, all the output .raw files with updated states from SUGAR-IMB based solver are added to the same repository.
The condensed results with comparison against standard non-linear optimization tools are shown in Table 1. The SUGAR-based approach with the IMB homotopy method demonstrates better convergence across all (15) scenarios compared to fmincon that converged for 12 out of 15 scenarios, IPOPT that converged for 11 out of 15 scenarios and MIPS that converged for 5 out of 15 scenarios. While the SUGAR-IMB homotopy solver does take more iterations compared to other tools to reach an optimal solution, the net runtime for the SUGAR-IMB solver is lower in many of those large scenarios when compared against other methods. This is primarily because the homotopy methods used in SUGAR-IMB ensure smooth convergence for each sub-problem. The standard commercial tools, however, perform more internal function evaluations between each iteration to converge to the final solution. Also, it's worth noting that SUGAR-IMB outperforms other methods for congested large networks (as these tend to have more complex solution space) as shown in the bottom two rows of the results table.

Conclusions
In this paper, we develop a methodology for robustly solving hard-to-solve AC-OPF problems especially when good initial conditions are not available. The proposed framework introduces a homotopy method called Incremental Model Building approach that initially modifies the problem to solve a trivial problem and then gradually constrains the problem to solve the AC-OPF exactly in the final homotopy factor. The homotopy method has its roots in circuit-theory and with inclusion of grid physics can satisfy critical criteria for convergence of homotopy methods. To demonstrate the robustness of the proposed framework, we show results for hard-to-solve networks and compare those against other standard non-convex optimization engines. The Incremental Model Building framework acts as a basis for future planning optimizations of the grid, such as optimal power flow and transmission expansion, that require a robust engine capable of scaling to large systems regardless of initial conditions.

Acknowledgement
This work was supported in part by the National Science Foundation under grant ECCS-1800812.