A novel MILP-Newton approach for power ﬂow analysis

This paper presents a novel mixed integer linear programming-Newton approach for solving constrained AC power ﬂow problems, based on a recently proposed linear programing-Newton method. The linear programming-Newton method overcomes some deﬁciencies of the classical Newton’s method, while retaining its quadratic convergence in the neighbourhood of a solution. The discrete nature of some control variables, such as transformer tap ratios and bus shunt susceptances, are often ignored in power ﬂow algorithms due to difﬁculties that they introduce. The violation of operational limits in power ﬂow algorithms is typically corrected in ex-post procedures, such as PV-PQ bus switching, which may lead to infeasibilities. Finally, under some circumstances, power ﬂow problems may not have a solution; in these cases, a desirable property of a power ﬂow algorithm would be the ability to obtain a nearly feasible solution. The approach proposed in this paper handles these difﬁculties in a uniﬁed framework and is able to deal directly with discrete variables and operational limits. When there are no power ﬂow solutions, a nearly feasible solution that minimizes a merit function is obtained. Numerical experiments show that the proposed approach is robust and obtains accurate solutions in the presence of discreteness conditions and operational limits.


INTRODUCTION
In order to meet the increasing electricity demand, power systems have grown both in size and complexity and, as a result, computational tools became mandatory for their safe and reliable operation. Among such tools, AC-Power Flow (AC-PF) algorithms are widely employed in power systems analysis. The arising of new technologies such as distributed generation and hybrid systems has presented new challenges for this research topic, making the power flow an important research topic even today [1,2]. Newton-like algorithms for AC-PF computations became standard in practice, because they often provide, fast and accurately, the steady state variables for the operator [3]. An example of such algorithms is the AC-PF toolbox available in the MATPOWER simulation package [4], which employs classical Newton's method and a polar form of power balance equations.
computations include use of techniques based on Runge-Kutta formulas [12,13], Levenberg's and Lyapunov-based methods [14], reformulations through mixed complementarity and optimization (interior point method) [15], as well as the factorized load flow of [16], which is an extension of the factorized solution methodology in state estimation to compute AC-PF solutions. Methods based on complex analysis are also available, such as the holomorphic embedding method [17], which employs the truncated Maclaurin series and Padé approximants to compute the power flows.
In AC-PF computations, it is necessary to take into account that some control variables, such as reactive power injections in generation buses, voltage magnitudes, and on-load tap-changer (OLTC) in-phase transformer tap ratios, need to be kept within their respective operational limits. Since classical Newton's method assume that variables are not bounded, power flow tools usually comprise two nested loops: an inner loop, where Newton's method is employed to reduce the power mismatches, and an outer loop to check control variables and enforce limits. However, this algorithmic scheme may slow down the convergence process or even diverge. For this reason, many AC-PF algorithms either completely ignore such limits or consider them only partially.
For example, in the AC-PF algorithm of MATPOWER, generator reactive power limits are enforced by fixing reactive injections at the limit for those generators with a violated reactive power limit, converting the corresponding bus into a PQ bus, and solving the AC-PF again, until no more violations are found. Nonetheless, this reactive power control procedure is performed at the expense of voltage setpoint, which may lead to inadequate voltage magnitudes.
Another issue that may arise when solving AC-PF problems is the fact that, in some situations, there is no solution within the operational limits. Heavy-loaded systems and occurrence of contingencies, such as a branch or generator outages, are likely to lead to unsolvable AC-PF problems. When this happens, a nearly feasible solution can be of great value for the power system operator. Several approaches described in the literature disregard this aspect and most Newton-like algorithms diverge when no solution is available. However, there is a growing concern about this aspect, and some papers such as [15] and [18] try to minimize power mismatches violation when a feasible solution cannot be attained.
Besides the aforementioned difficulties, most AC-PF solution approaches also disregard that some control variables, such as OLTC in-phase transformers and bus shunt susceptances can only be adjusted to a predefined set of discrete values. These variables are typically fixed or assumed to be continuous when their adjustment is needed.
In order to recover a practical feasible solution, one can round-off discrete variables to the nearest discrete feasible value, and run the PF algorithm again. However, there is no guarantee that a feasible solution will be obtained after this rounding-off procedure. Although the importance of discrete variables has gained the attention of researchers in the context of optimal power flow (OPF) studies [19,20], they are still neglected in AC-PF solution approaches.
In this paper, we propose a novel approach to solve constrained AC-PF problems, which is based on the linear programming (LP)-Newton algorithm proposed by [21] and extended by [22]. When compared to classical Newton's method, the LP-Newton method features several advantages: it is able to solve systems of equations with n ≠ m, Ω ≠ ℝ n , complementarity equations and piecewise smooth equations. It also converges quadratically to a solution of (5) under weak assumptions [21]. The main step of the LP-Newton method requires the solution of a LP problem to compute the search directions to update the variables.
In order to handle constrained systems of equations with discrete variables, we propose a mixed-integer LP (MILP)-Newton approach. In our approach, LP subproblems that lie at the core of LP-Newton method are transformed into MILP problems, so that search directions preserve feasibility of discrete variables. Advantages and disadvantages of this approach are discussed, as well as some practical details that may improve its efficiency.
We assess the performance of our approach by performing numerical tests with IEEE power systems of 14, 30, 57, 118 and 300 buses. Numerical results indicate that our approach is able to obtain AC-PF solutions when they exist or minimize mismatch violation if a solution is not available while keeping continuous control variables within the specified limits and assigning adequate discrete values to discrete variables. The main contributions of this paper can be summarized as follows: i. we present a globally and quadratically convergent LP-Newton method, which is robust and able to compute steady-state variables in AC-PF while respecting operational limits, without the need of iterative PV-PQ bus switching. It is also able to obtain minimum power mismatch solutions when a solution within operational limits possibly does not exist; ii. we propose a MILP-Newton approach, which is an extension of the LP-Newton method that is able to obtain solutions for AC-PF problems with discrete control variables, by solving MILP subproblems; iii. we devise some modifications and heuristics that may improve the practical performance of the MILP-Newton approach.
In addition to the main contributions presented above, Table 1 summarises other differences between the presented paper and other approaches proposed in the literature. The asterisk in this table is due to the fact that, while our paper does not present numerical results for a constrained system of equations with complementarity constraints, such constraints can be handled by the LP-Newton method [21].
In addition to this introduction section, this paper is organized as follows: in Section 2 we present the mathematical model of a constrained AC-PF problem. Section 3 presents the basic local and global convergence theory of LP-Newton method, as well as a novel MILP-Newton approach for solving constrained systems of equations with discrete variables.

CONSTRAINED AC POWER FLOW
In a typical AC-PF computation, one seeks to calculate the steady-state operation voltage magnitudes and angles. This is done in a two-step procedure: first, classical Newton's method is applied to solve the nonlinear system of equations defined by active power balance equations of PV and PQ buses, and reactive power balance of PQ buses, with a fixed slack bus angle. Once a solution is found, it is used to compute the remaining unknown quantities, such as reactive power generation. However, this procedure has some caveats. Firstly, Newton's method may diverge, especially in ill-conditioned cases. Even if convergence is attained, it is not possible to ensure that the variables will be kept within the required operational limits.
Based on these facts, in this paper, the AC-PF problem is formulated as a constrained system of equations as follows: Active and reactive power flows from bus k to bus m are functions of the variables v, and t , given by the following expressions [23]: In the AC-PF model (1), active power balance equations for PV and PQ buses (1a) and reactive power balance equations for PQ buses (1b) are the ones typically employed in Newton-like algorithms for AC-PF computations [3,18,20,23,24]. Reactive power production by PV buses and operational limits (voltage (1d) and reactive power production limits (1e)) are often enforced in ex-post procedures [3,23,24], which may cause infeasibilities that cannot be corrected by running the Power Flow (PF) algorithm again. To prevent such issues, the development of algorithms that can deal with the constrained system of equations (1) as a whole is important.
A continuous constrained AC-PF model is obtained if the discretization constraints (1f) and (1g) are relaxed as follows: For the sake of simplicity and to meet the notation of next section, the constrained system of Equation (1) can be restated in the compact general form by setting and where: = (v, , q, t, b sh ) ⊆ ℝ n is the vector of variables; n = || + | ′ | + || + | | + | sh | is the number of variables; v and v are the lower and upper bound vectors for voltage magnitudes, v; q and q are the lower and upper bound vectors for reactive power generation, q; ∏ (k,m)∈  t km is the Cartesian product of the discrete sets for in-phase transformer tap ratios, t ; and ∏ k∈ sh  sh k is the Cartesian product of the discrete sets for bus shunt susceptances, b sh .

PROPOSED METHOD
The numerical solution of constrained systems of equations of the form (5), where Ω ⊆ ℝ n is a non-empty closed set and F : ℝ n → ℝ m is a continuous function, is of paramount importance in many applications. Examples include AC-PF computation and Karush-Kuhn-Tucker (KKT) systems arising from optimization problems, such as the OPF problem.
The most known case of the constrained system of Equation (5) happens when n = m and Ω = ℝ n , which corresponds to solve a square system of equations. Algorithms for this situation are often based on classical Newton's method, due to its quadratic convergence when the initial iterate is within a proper neighbourhood of a solution. Despite its numerical efficiency, Newton's method may diverge if the initial iterate is not chosen properly or the problem has no solutions. Additionally, classical Newton's method is not well-suited to handle limits on the variables, .

LP-Newton method
The LP-Newton method was initially proposed in [21] to solve the constrained system of Equation (5). We will assume that the solution set is non-empty: Given an iterate k ∈ Ω, the search direction employed to update the variables, d k , is calculated by solving the following subproblem [21] in the variables (d, ) ∈ ℝ n × ℝ: where G ( k ) denotes the Jacobian matrix of F at k (or a proper generalized Jacobian, if F is non-smooth). From a theoretical point of view, one can use any norm in (9) and any closed set Ω. However, it is important to notice that suitable choices of norm and set Ω can make the solution of this subproblem significantly easier. For example, if Ω is a convex set, then (9) is a convex optimization subproblem. One should notice that subproblem's (9) goal is to minimize the additional variable . Hence, whenever we say that d k is a solution of (9), it is implied that d k , alongside a k = ( k ), solves the subproblem. The following result shows that (9) has a solution for any choice of k ∈ ℝ n . Proposition 1 ([21]). For any k ∈ ℝ n , the subproblem (9) has a solution, and the optimal value of (9) is zero if and only if k is a solution of (5).
Due to Proposition 1, the following basic LP-Newton algorithm is well-defined for any choice of starting point 0 ∈ Ω. Algorithm 1 is not only well-defined. As shown in [21], it also has very strong local convergence properties and converges quadratically under assumptions that are weaker than some standard conditions proposed in the literature to establish local convergence of other Newton-like methods.

Linesearch globalization
Besides the fast local convergence, robustness is another important feature for an iterative algorithm. Thus, global convergence, that is the ability to converge regardless the choice of initial point, is a desirable property for algorithms. In most cases, global convergence is established using either a linesearch or trust-region approach, and a merit function is employed to gauge the progress that a direction provides towards a solution of problem.
In the case of unconstrained optimization, the obvious choice of merit function is the objective function itself. For non-linear systems of equations, several choices of merit functions are available, each with some advantages and disadvantages [25]. Ideally, characteristics of subproblems that define search directions should hold an intimate connection with the chosen merit function [26].
To present the globalization scheme described in [22], we will assume throughout this subsection that ‖ ⋅ ‖ denotes the ∞norm and that Ω is a non-empty polyhedral set (which is the case for continuous relaxations of the AC-PF presented in Section 2). Under these assumptions, subproblem (9) is just a LP problem, whose solution can be efficiently obtained by several solvers. In this case, the natural choice of merit function is given by, which is a violation measure for the constrained system of equations (5), and inherits the same norm used to solve the ALGORITHM 2 Globally Convergent LP-Newton Method 1. Choose ∈ (0, 1), ∈ (0, 1), 0 ∈ Ω, > 0 and set k = 0. 2. If m( k ) ⩽ , then STOP. Return k . 3. Solve (9) to obtain (d k , k ) and set = 1. 4. While (13) does not hold, set = . 5. Set k = , update the variables by k+1 = k + k d k , increase the iteration counter k by 1 and return to Step 2. subproblems. It is clear that if (5) has a solution, then it is also a solution for the minimization problem: whose necessary optimality condition is given by, where m( ⋆ ) is the Clarke generalized gradient of a Lipschitz continuous function m near ⋆ and N Ω ( ⋆ ) is the normal cone to Ω at ⋆ [27]. However, the converse is not true: some local minimizers of (11) satisfying the necessary optimality condition (12) may not be solutions for (5).
Once the search direction d k is computed as the solution of (9), a backtracking linesearch is performed to ensure that a sufficient reduction of the merit function is achieved, by checking the Armijo condition [25], given by where ∈ (0, 1] is the trial step-length, ∈ (0, 1) is a parameter and is an estimate of the directional derivative of m at k along the direction d k [22]. Therefore, a globalized LP-Newton algorithm can be obtained by simply adding a linesearch procedure with a merit function after the computation of search directions. Algorithm 2 summarizes these steps.
The following theorem, proven in [22], establishes the global convergence of Algorithm 2. Theorem 1 ([22]). Algorithm 2 is well-defined and, for any starting point 0 ∈ Ω, it either terminates finitely with an iterate k satisfying 0 ∈ m( k ) + N Ω ( k ) or generates an infinite sequence { k } and any accumulation point ⋆ of this sequence satisfies (12).
Remark 1. Theorem 1 ensures global convergence of Algorithm 2 in the sense of finding a point satisfying the necessary conditions (12). If such point does not satisfy m( ) = 0, then it is not a solution for the constrained system of equations (5). Despite this fact, such result is still highly relevant, since it ensures that if the sequence of iterates have an accumulation point, then it locally minimizes the merit function, even when there are no solutions for (5).

Handling discrete variables: MILP-Newton approach
Proper handling of discrete variables is a challenging task for any algorithm designed to solve (5), because convexity assumptions cannot be made. This means that the globally convergent Algorithm 2 cannot be directly applied in this context. Nevertheless, local quadratic convergence is still preserved, as long as the required assumptions are satisfied and Ω is a closed set.
In order to solve constrained systems of equations with discrete variables, we propose to exploit the structure of each LP-Newton subproblem, by incorporating discreteness conditions as constraints. The search directions which are calculated ensure that iterates remain within limits in the case of continuous variables and are discrete in the case of discrete variables.
Let  be the set of indexes of the continuous variables in ,  be the set of indexes of discrete variables in , which can assume values in an evenly spaced discrete set, and  be the set of indexes of discrete variables in , which can assume values in an unevenly spaced discrete set. Therefore, we propose to solve the following subproblem to compute search directions: where  i is the non-empty finite set of discrete settings for a discrete variable indexed by i ∈  ∪  . Here, for the sake of clarity, we present the discreteness condition separately in (15e). Such constraints are reformulated as algebraic equations, by introducing extra integer variables, resulting in a MILP problem.
It is important to notice that there is no significant difference between these subproblems, since (15e) can be regarded as a part of (15d). We only stated the subproblem in such manner to emphasize that discrete variables need appropriate techniques to be handled. Therefore, the Jacobian matrix G is computed in the same way for both the continuous and discrete AC-PF and is similar to those employed in classical AC-PF tools. Its columns are the partial derivatives of the components of F with respect to v, , q, t and b sh .
One should be aware that using the bus conductance (G ) and bus susceptance (B) matrices to calculate these derivatives is not advisable in the model presented in this paper, because transformer taps are variables. As such, these derivatives were calculated based on Equations (2) and (3). What separates the LP-Newton method from the MILP-Newton approach is the fact that the latter one requires additional constraints to enforce that discrete variables are calculated at each step.

Evenly spaced discrete variables
Let i ∈ . We can assume that | i | ⩾ 2 because, otherwise,  i is a singleton and i can be regarded as a constant instead of a variable. In this case, the step h i between two consecutive elements in  i can be calculated as Therefore, for i ∈ , constraint (15e) can be stated as: where u i is an extra integer variable, added for each i ∈ .

Unevenly spaced discrete variables
In this case, we can state constraint (15e) as: where u i, j are extra binary variables, that are assigned to each possible discrete setting of each unevenly spaced discrete variable. Equation (18b) is a logical constraint which ensures that a single discrete value is selected in each set  i .
The main reason why we cannot employ Algorithm 2 in this case is because, if the step is shortened during the linesearch FIGURE 1 LP-Newton and MILP-Newton algorithm flowchart: the dashed arrow indicates that the MILP-Newton approach must skip the linesearch procedure, then the discreteness condition is not preserved. As a result, when solving problems with discrete variables we only perform full-steps, as in Algorithm 1. For the sake of clarity, Figure 1 presents a flowchart for these algorithms. The dashed arrow indicates that the MILP-Newton approach skips the linesearch to ensure that discrete variables are obtained at each iteration.

3.4
Practical details

Advantages and disadvantages
The main disadvantage of the LP/MILP-Newton approaches is the need to solve LP/MILP subproblems. It is obvious that the cost per iteration is higher than traditional Newton's method, which only relies on the solution of linear systems of equations. However, this increased cost is paid off by the possibility to solve many problems that Newton's method cannot. An example is the case of problems with complementarity equations. If xy = 0 is an equation of the system, then the corresponding Newton equation is y k Δx + x k Δy = −x k y k . Thus, if an iterate becomes zero, say, x k = 0 and y k ≠ 0, the equation will become y k Δx = 0 which, in turn, implies that Δx = 0. As a result, the variable x is no longer updated, and the method gets stuck. On the other hand, if both x k = 0 and y k = 0, then the Jacobian matrix is not of full rank, which also lead to numerical issues. Anyway, Newton's method is likely to fail.
In contrast, LP/MILP-Newton approaches have the innate ability to: solve problems with additional inequalities (defined by the set Ω), complementarity equations and piecewise smooth systems equations with non-isolated solutions.

Constrained AC-PF subproblem
Under the assumption that ∞ -norm is used, the LP-Newton method subproblem (9) can be restated as: where I n is the identity matrix of order n and e represents vectors of ones with appropriate dimensions. If Ω is polyhedral, then this is a LP problem, which can be efficiently solved.
In the case of the continuous relaxation of the AC-PF problem presented in this paper, (20f) represents the upper and lower bound constraints on variables, . Considering that = (v, , q, t , b sh ) T and = (v, , q, t , b sh ) T are, respectively, the upper and lower bound vectors for the variables, this constraint can be rewritten as, Assuming that the search direction is partitioned according to , that is, As such, constraint (21) simply provides upper and lower bounds on the search direction components, which can be easily set in LP solvers.
In the case with discrete variables, we simply add extra integer or binary variables and equality constraints to this subproblem, according to Section 3.3. For in-phase transformer tap ratios, which are evenly spaced discrete variables, we introduce the constraints: where h t i j is the step between two consecutive tap ratios. For bus shunt susceptances, which are unevenly spaced discrete variables, we introduce the constraints:

Non-monotone linesearch
For the globally convergent Algorithm 2, practical performance can often be improved by employing a non-monotone variant [28]. This is done replacing the Armijo condition (13) by: where is a positive integer parameter that indicates the number of reference values of the merit function that are stored for comparison. Clearly, if = 1, the usual monotone Armijo linesearch condition (13) is recovered. As increases, the linesearch becomes less conservative, and is more likely to accept a fullstep, even if it increases the merit function value.

Subproblem truncated solutions
LP subproblems are usually solved very fast with current stateof-art solvers. However, for MILP subproblems, a zero gap solution may take a long time to be obtained. In these situations, it is convenient to interrupt the solver when a feasible solution is found (the incumbent solution) and a fixed time threshold is reached. This can be done, as shown in the general Newton framework developed in [29], which includes the LP-Newton method. Even with truncated feasible solutions, local quadratic convergence is still attained under the required assumptions. This means that, in some situations, a subproblem does not need to be solved until optimality. However, for the globally convergent Algorithm 2, optimal solutions for the subproblem are required to compute an estimate for the directional derivative.

Bad scaling
It may happen that, due to bad scaling, the LP solver used for (9) fails to converge. In this situation, one can perform the variable change = ‖F ( k )‖, and try to solve the subproblem again. Although this is not the case for numerical experiments presented in this paper, adding this procedure may increase robustness of the algorithm. Further implementation details can be found in [22].

AC-PF computation
If a discrete variable has a high number of possible discrete states or there are many discrete variables, the performance of MILP solvers can be degraded. This is the case of OLTC inphase transformer tap ratios in this paper (see Section 4). To improve the MILP solver performance in these cases, we compute search directions in such a way that a transformer tap variable can either retain its current value or be updated to the discrete setting that is immediately above or below its current value.

NUMERICAL RESULTS
In order to evaluate the performance of LP/MILP-Newton approaches to solve constrained AC-PF problems, tests were performed with IEEE-14, 30, 57, 118 and 300 buses power systems from MATPOWER package. The algorithm was implemented in MATLAB 2012a and LP/MILP subproblems were solved using MATLAB's toolbox from CPLEX 12.6.0. In order to obtain accurate search directions, the solver stopping accuracy was set to 10 −9 .
Tests were performed on a machine with an Intel Core i9-9900 @ 3.60GHz processor and 16GB of RAM. Lower and upper bounds on voltage magnitudes were obtained from MAT-POWER test cases , i.e. 0.94 and 1.06 p.u. to the lower and upper limit values, respectively. The discrete set for transformer tap ratios,  t km , ranges from 0.88 to 1.12 p.u. with evenly spaced discrete steps of 0.0075 p.u., totalizing 33 possible discrete tap settings, for all transformers. The discrete sets for bus shunt elements,  sh k , are presented in Table 2. The tolerance for convergence was set to = 10 −6 , while the parameters of linesearch in Algorithm 2 were = 0.001 and = 0.5. Regarding the starting point, flat start was used for voltages, the average value of upper and lower limits for reactive power generations, while transformer tap ratios and bus shunt susceptances were started using the test case data.

Continuous relaxation and linesearch parameter influence
We first analyse how influential is the non-monotone linesearch parameter . To do so, we solve continuous relaxations of all test cases, by considering (4), with several choices of . Table 3 summarizes the results obtained using Algorithm 1, which does not perform a linesearch, and the globally convergent Algorithm 2, using different choices for the non-monotone linesearch parameter ( = 1 corresponds to a monotone linesearch). Figure 2 shows the convergence behaviour in numerical experiments with the IEEE-300 bus case and Figure 3 shows the voltage magnitudes at the solutions.
These results indicate that the non-monotone linesearch usually outperforms the pure LP-Newton method (Algorithm 1) and its globally convergent monotone version (Algorithm 2 with = 1). This is related to the fact that the non-monotone linesearch allows for more full LP-Newton steps to be taken, which accelerates the convergence, while avoiding significant oscillations of the merit function, as in the pure LP-Newton algorithm. On the other hand, the monotone Armijo condition is more stringent, forcing the algorithm to shorten the step-length more often.
For the IEEE-300 case, convergence trajectory is the same for = 2 and = 3. It is important to notice that more freedom to the non-monotone linesearch, i.e., higher values of , does not necessarily translate into a reduction of iterations or computation time. Nevertheless, all the algorithms are able to obtain a solution, within the operational limits, which is a feature of the proposed algorithms in this paper.

Robustness and importance of oltc in-phase transformer tap ratios and bus shunt susceptances as variables
In this subsection we show that considering transformer tap ratios and bus shunt susceptances is very important to obtain constrained AC-PF solutions with the required accuracy. We consider three cases for the IEEE-300 buses power system:  Table 4. As shown in Table 4, when transformer tap ratios and bus shunt susceptances are fixed, it is not possible to find a solution with a mismatch below the required tolerance, and the algorithm stops when it reaches a maximum number of iterations of 300. Despite this fact, the solution in this case may still be useful for power system operators because the maximum mismatch is of 1.3366 × 10 −3 p.u. and all variables are within the required operational limits. This particular test also demonstrates that our approach is robust and able to obtain solutions that minimize ‖F ( )‖, even when a solution possibly does not exist.
In the discrete case of this experiment, a higher computation time is observed due to the need to solve MILP subproblems. Nevertheless, an adequate solution with discrete transformer tap ratios and bus shunt susceptances is found. Figures 4 and  5 show the convergence behaviour and voltage magnitudes at solutions in these tests.
Regarding Figure 5, it is possible to notice a higher voltage magnitude profile in the fixed case. This naturally suggests that infeasibilities in this situation happen because it is impossible to keep voltages within its operational limits and satisfy power balance equations at the same time. Indeed, if we slightly increase maximum voltage magnitudes to 1.07 p.u., our approach is able to converge in 12 iterations in the fixed case. On the other hand,

LP-Newton × classical Newton's AC-PF algorithm
In order to put our approach into perspective, we also compare it to Newton's method implemented in MATPOWER's AC-PF routine. The goal of such comparison is to show that both approaches have their strengths and weaknesses. For the purpose of comparison, we enforce generators reactive power limits in MATPOWER's algorithm.
Results in Table 5 indicate that, as expected, classical Newton's method is usually faster than our approach, because it only solves linear systems of equations to compute search directions. However, it lacks the desirable ability to keep variables within their operational limits. For example, in the IEEE-300 test case, lower and upper voltage magnitude limits are violated at some buses, as shown in Figure 6. Additionally, with the proposed MILP-Newton approach, discrete solutions within such limits can also be calculated. As such, our approach is advantageous when these constraints are relevant for the power system operator.

MILP-Newton starting point
Finding discrete solutions in AC-PF problems can be a difficult and time consuming task. As shown in Table 5, for the IEEE-300 bus test case, a solution with discrete variables was found in 96.4825s. To improve computation time, it is possible to use the solution of the continuous relaxation obtained by the LP-Newton method as a starting point for the MILP-Newton approach. However, there is one detail that needs to be taken into account. The local convergence theory and algorithms for LP/MILP-Newton approaches require that a starting point must satisfy 0 ∈ Ω. This means that the MILP-Newton algorithm should, ideally, be started with discrete values for discrete variables. For obvious reasons, the continuous relaxation solution of LP-Newton method is unlikely to provide such starting point.
A possible way to handle this situation is to round-off discrete variables to the nearest allowed discrete value, an approach that we call LP+RO+MILP-Newton. However, we found out in our experiments that using the continuous solution itself as starting point can also provide good results in terms of time savings and iterations. Since the MILP-Newton approach preserves feasibility of the iterates in relation to Ω, after the first iteration, all   discrete variables will immediately assume discrete values, and feasibility will be maintained throughout the remaining iterations. This starting strategy is called LP+MILP-Newton.
The results of such strategies are presented in Table 6. These results demonstrate that significant time savings can be obtained by starting the MILP-Newton approach with good quality solutions, due to the reduction of MILP subproblems that are solved. It also should be emphasized that there is not a clear winner between these strategies: while LP+RO+MILP-Newton had better results for the IEEE-118 test case, LP+MILP-Newton performed better for the IEEE-57 and IEEE-300 test cases. Figure 7 shows that all transformer tap ratios are properly discretized in the IEEE-300 bus test case, as well as the continuous relaxation transformer tap ratios obtained by Algorithm 2 with = 2. In Figure 7, dashed black lines are the allowed discrete tap ratios. Table 7 shows that bus shunt susceptances are also properly discretized. It should be noticed that the MILP-Newton, LP+RO+MILP-Newton and LP-MILP-Newton approaches do not simply round-off the results obtained by the continuous relaxation.

Standard voltage magnitude limits
The results presented in previous sections consider the upper and lower bounds on voltage magnitudes of 1.06 p.u. and 0.94  Table 8. For this set of tests, the LP+MILP-Newton had a better performance than the LP+RO+MILP-Newton, since the latter did not converge to a discrete solution of the IEEE-300 case. This indicates that using the rounded-off continuous solution as starting point may degrade the performance of the proposed MILP-Newton approach. Despite this fact, discrete solutions were obtained with both MILP-Newton and LP+MILP-Newton.
The voltage magnitude profile for the IEEE-300 test case solution is presented in Figure 8. Transformer tap ratios and  bus shunt susceptances for both, continuous and discrete simulations, are presented in Figure 9 and Table 9, respectively.

Comparison with previous works
As far as we are aware, there are no recent papers that take into account the same operational constraints and discrete variables that we consider in the AC-PF model of our work. We chose two  papers that also consider voltage magnitude limits and reactive power generation limits for comparison. However, this comparison is only to put our approach in perspective, since there are some differences between the models. The paper [15] presents results for a continuous AC-PF problem, with complementarity constraints. The authors do not consider in-phase transformer tap ratios and bus shunt susceptances as variables. Among the continuous optimization solvers employed in that paper, the interior-point solver IPOPT was the one with overall better performance. For this reason, we compare it to the results of the globally convergent (continuous) LP-Newton method presented in Table 8. This comparison is available for all the test cases, in Table 10. The numerical experiments show that the LP-Newton method finds a solution within a reasonable amount of time.
In [17], the authors propose the Holomorphic Embedding Method (HEM) for power flows. In this algorithm, after performing bus-type switching checks, the voltage  [17] also indicate that is necessary to investigate more elegant and efficient methods for bus type switching and tap changing that take advantage of the HEM. Therefore, we compare the LP+MILP-Newton approach proposed in our paper with [17], by computing the execution time ratio of the LP+MILP-Newton method presented in Table 8 and Classical Newton's method implemented in MATPOWER presented in Table 5. This comparison is available for the IEEE-118 and 300 buses test cases, in Table 11.
These results show that the LP+MILP-Newton approach had a higher execution time ratio. This was expected due to bus shunt susceptances, which are significantly harder to discretize in comparison to transformer tap ratios. This fact was noticed long ago, for example, by [30], in the context of Optimal Power Flow problems.
According to these authors, the most important discrete variables are for shunt capacitors and reactors, but transformers are also of some importance. Rounding-off to the nearest step is marginally acceptable for controls such as transformers whose steps are small and uniform in size. But the errors of rounding are quite large for controls whose steps are large and nonuniform.
This is certainly the case of our paper: while transformers have discrete step sizes of 0.0075 p.u., bus shunt suscepances can have a discrete step of 4.5 p.u., according to Table 2. As a result, if there is a large bus shunt susceptance change between two iterations, we can expect a noticeable change in some mismatches.

LP-Newton local infeasibility detection
In this set of simulations, we consider the IEEE-118 buses test case to examine the behaviour of LP-Newton method in situations that a local solution for the AC-PF, within the operational limits, possibly does not exist. To do so, we consider voltage magnitude limits as ±5%, ±4%, ±3%, ±2% and ±1% of the rated voltage (1 p.u.). We only consider the continuous case in this simulation for two reasons: the continuous case algorithm should, ideally, be executed before any attempt of solving the discrete case, because a discrete solution only exists if there is also a continuous solution; secondly, the theory of the globally convergent LP-Newton method ensures that if an iterate k is such that F ( k ) ≠ 0, i.e., if the mismatches are not zero, then Δ( k ) = 0 if, and only if, 0 ∈ m( k ) + N Ω ( k ), where m( ) = ‖F ( )‖ (see Lemma 3.2 and the discussion below its proof in [22]).
As a consequence, if we identify that the LP-Newton method does not stop by its usual criterion, ‖F ( k )‖ ⩽ , but |Δ( k )| is sufficiently small, we can conclude that the method is converging to a point that locally minimizes the merit function, but is not a solution for the problem. As such, to perform these experiments, we included the additional stopping criterion |Δ( k )| ⩽ 10 −4 . We adopted a higher tolerance here, otherwise it would take several iterations to detect local infeasibility. The results are summarised in Table 12.
These results show that, when the method is not able to find a solution for the continuous AC-PF, the value of |Δ( )| gives a valuable information about local infeasibility. For example, in the case of voltage magnitudes ranging from 0.99 p.u. to 1.01 p.u., the maximum mismatch was of 0.1053 p.u., with |Δ( )| = 9.8372 × 10 −5 , which suggests that there are no solutions for the problem. Figure 10 shows that, in this situation, the maximum mismatch reach a steady value around iteration 50, while the value of |Δ( )| decreases until infeasibility is detected. As we widen the voltage magnitude interval, the method finds a solution in the interval ranging from 0.97 p.u. to 1.03 p.u.

CONCLUSION
This paper proposed a novel MILP-Newton approach for solving constrained systems of equations with discrete variables, based on the recently proposed LP-Newton method [21]. It is able to solve systems of equations that a classical Newton's method either fails or is not well-suited to handle. This is the case of problems with complementarity equations, inequality constraints, piecewise smooth functions, non-square systems or systems with discrete variables. Such systems of equations happen naturally when solving AC-PF problems or formulating the KKT conditions of an optimization problem. Each iteration of this approach requires the solution of LP subproblems (in the continuous case) or MILP subproblems (if there are discrete variables) to compute search directions, as long as the infinity norm and a polyhedral constraint set are employed. Numerical experiments on constrained AC-PF problems with discrete control variables, such as OLTC in-phase transformer tap ratios and bus shunt susceptances, showed that the method obtains solutions of high accuracy in a good computation time within the required operational limits, which is a lacking feature of most AC-PF algorithms. We believe that our approach holds potential to solve more representative AC-PF models. In future work, we plan to extend its applications to AC-PF models with complementarity constraints as well as optimization problems.