{'fraction_non_alphanumeric': 0.028452227908894658, 'fraction_numerical': 0.007131248938802241, 'mean_word_length': 4.732480533926585, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reduced-Space Interior Point Methods in Power Grid Problems

Due to critical environmental issues, the power systems have to accommodate a significant level of penetration of renewable generation which requires smart approaches to the power grid control. Associated optimal control problems are large-scale nonlinear optimization problems with up to hundreds of millions of variables and constraints. The interior point methods become computationally intractable, mainly due to the solution of large linear systems. This document addresses the computational bottlenecks of the interior point method during the solution of the security constrained optimal power flow problems by applying reduced space quasi-Newton IPM, which could utilize high-performance computers due to the inherent parallelism in the adjoint method. Reduced space IPM approach and the adjoint method is a novel approach when it comes to solving the (security constrained) optimal power flow problems. These were previously used in the PDE-constrained optimization. The presented methodology is suitable for high-performance architectures due to inherent parallelism in the adjoint method during the gradient evaluation, since the individual contingency scenarios are modeled by independent set of the constraints. Preliminary evaluation of the performance and convergence is performed to study the reduced space approach.

The real-world transmission lines are limited by the instantaneous amount of power that can flow through the lines due to the thermal limits. The apparent PFs of the transmission lines, S f ∈ N L and S t ∈ N L , are therefore limited by the power injections at both "from" and "to" ends of the line, which cannot exceed prescribed limits F max L . Squared values of the apparent PF are usually used in practice.

Optimal Power Flow
Since the formulation of optimal power flow (OPF) by Carpentier [9] as a continuous nonlinear programming (NLP) problem, OPF has become one of the most important and widely studied constrained nonlinear optimal control problems. It is concerned with the optimization of the operation of an electric power grid subject to physical constraints imposed by electrical laws and engineering limits. The objective is to identify the operating configuration that best meets a particular set of evaluation criteria. These criteria may include the cost of generation, transmission line losses, and various requirements concerning the system's security, or resilience with respect to disturbances. The standard OPF problem is formulated as minimization of the active power generation cost, subject to a set of equality and inequality constraints: The objective function is a quadratic functional of the active power generation, with the cost coefficients a , b , c defining the cost function of each generator. The AC nodal PF balance equations for each bus in the grid represent the nonlinear equality constraints (1b), g ∈ 2N B , where the exponential terms in the polar voltage representation are replaced using the Euler identity, and subsequently, the complex quantities are split into their real and imaginary components, g = (ℜ(g), ℑ(g)). Line power flows limits h ∈ 2N L are imposed by a set of inequality constraints (1c), where Over the last five decades almost every mathematical programming approach that can be applied to OPF has been attempted. Sequential linear programming (LP) has been presented by Stott and Hobson [32] and Alsac et al. [2]. Despite its performance it introduces inaccuracies and might lead to infeasible solutions. Sequential quadratic programming (SQP) was applied by Burchett [5] delivering accurate and feasible solutions but the computational times are relatively high, rendering the whole solution approach inappropriate for large-scale power systems. IPMs have became a successful tool for solving constrained optimization problems. The origin of the IPMs popularity reaches back to the 1984 when Karmarkar [22] announced a polynomial time linear program that was considerably faster than the most popular simplex method to date. Furthermore, IPMs can be applied also to quadratic and other nonlinear programs, unlike the simplex method which can be applied only to linear programming. The main advantages of the IPM lie in the ease of handling inequality constraints by logarithmic barrier functions, speed of convergence, and a fact that strictly feasible initial point is not required. Linear and nonlinear IPM were applied to the OPF problems since the 90s by Vargas et al. [34,18,33]. Several extensions of these algorithmic approaches are reported in recent papers; see, e.g., [16,14,24]. Although the LP-based approaches can already be applied close to real time for various practical OPF problems at least in medium size systems, further improvement in tractability need to be envisaged to get closer to the real-time AC OPF dream in large-scale systems, as pointed out by [8]. Another advantage of IPMs is that they are applicable to large problems and can easily exploit problem structure. The underlying linear algebra components of structured problems allow for large-scale parallel optimization problems on high-performance computers.
Modern trends in power grid operations and modelling render approximation based optimization techniques less attractive for coping with stressed operation conditions. One important advantage of NLP for OPF is that it naturally captures nonlinearities in stressed power system behavior, rendering them into an excellent tool for modelling and simulation of modern power system operations. As a result there is a great demand for new algorithms and software tools able to address strong nonlinearities in system behavior, in order to guarantee reliable and economic system operation. Additionally, in order to comply with current and near future power systems operating conditions, the incorporation of flexible security criteria and the consideration of fast control actions due to growing renewables penetration should be the natural, very ambitious and very challenging targets for the first generation of real-time AC OPF under uncertainty tools.

Security Constrained OPF
The security constrained OPF (SCOPF) [25], is an extension of the OPF problem, which finds an optimal operational state but at the same time takes into account a set of security constraints arising from the operation of the system under a set of postulated contingencies. It guarantees the whole power system can work under the nominal long-term cost-efficient operation plan, but can also remain in the operational state when some of the predetermined contingencies occur. However, each additional contingency corresponds to an additional set of constraints in the OPF problem specified for the associated power grid. The nominal scenario and all contingency states are coupled, rendering the whole problem computationally intractable for realistic size grids.
The SCOPF optimization problem considered in this work is the "preventive" SCOPF, although the same algorithmic improvements would apply also to the "corrective" variant. SCOPF problem is formulated as: Note that the SCOPF problem replicates the OPF constrains and variables for each contingency scenario c . The values of the non-automatic control variables are the same in all system scenarios, as expressed by the two non-anticipatory constraints (3f) and (3g). These declare that the voltage magnitude and the active power generation at the PV buses P V (also known as generator buses, the active power and voltage magnitude are specified, therefore the name PV bus) should remain the same as in the nominal scenario c 0 , regardless which contingency they are associated with. The only generator that is allowed to change its active power output is the generator at the reference bus r e f (also known as slack or swing bus) as its active power generation can be modified to refill the power transmission losses occurring in each contingency c . This implies that part of the optimization vector x will be shared between the scenarios and part of it will be local to each contingency. Therefore, the vector of variables can be partitioned into local components x c for each contingency ∀c ∈ {c 0 , c 1 , . . . , c N c } and the global (shared) part x g : This ordering of the variables allow implementation of the efficient structure exploiting algorithms described in the following sections. An important question arises in the SCOPF problem studies, which is how to select a reduced number of equivalent credible umbrella contingencies and associated variables that the SCOPF problem has to consider in order to obtain the same or nearly the same solution as with the full set of contingencies. In other words, how to identify contingencies that do not pose any security risks and at the same time do not restrict the optimal solution, therefore do not have to be considered in the SCOPF problem. The set of umbrella contingencies for a given SCOPF is strongly dependent on the operating conditions of the power grid (loading, grid configuration, etc.). Hence, as the parameters of the problem change, the membership in the set of umbrella contingencies also varies. Machine learning techniques are usually used for this purpose, including three classes of learning methods, namely machine learning, artificial neural networks and statistical pattern recognition [13,4]. However, the resulting set might be still very large and sophisticated optimization algorithms are needed. For the purpose of this study, such contingencies are selected that result in non-empty feasible region of the SCOPF problem (3). This restricts the selection of the transmission lines that can experience failure and do not compromise operation of the power grid. Such lines are characterized by the following properties: (i) no islands and isolated buses appear in the grid after the line failure, (ii) the reduced grid remains feasible in the PF sense, and (iii) only limited reactive power generation violations are allowed after the contingency occurrence. These considerations are usually made in the planning stage of the transmission grid design, but since synthetic power grid networks are used in this study, these considerations might not be satisfied for all conceivable contingencies.
Several techniques have been proposed in the literature aiming to reduce the computational complexity of the SCOPF problem, such as direct current (DC) approximations [15], filtering of contingencies using prior knowledge about the grid [6], [7] or nonlinear decomposition of the grid into zones with the aid of additional user information [20]. Other algorithms adopt the dual decomposition or the alternating direction method of multipliers, e.g., applied to the DC OPF problems in [10], [21]. Recently, in addition to classic CPU-based computing, graphical processing units (GPUs) have also been explored for security analysis in [11] or [36]. Alternatively, the SCOPF problem can be decomposed on the linear level. Interior point methods (IPMs) have been the most robust and successful tools for large-scale nonconvex optimization, and IPMs can easily exploit the problem structure [17]. Most of the computational time in IPMs is spent in the solution of the linear system, arising from the linearization of the Karush-Kuhn-Tucker (KKT) optimality conditions, the size of which grows with the number of considered contingencies and the size of the power grid. The KKT system is usually very large but sparse. Direct sparse linear solvers are the standard choice employed for the solution of the KKT. The solution can be accelerated using techniques tailored to the particular structure of the associated KKT linear system.
The parallel distributed Schur complement framework was applied to AC SCOPF problems and associated KKT systems in [19]. The work focuses on solving the Schur complement equations implicitly using a quasi-Newton preconditioned conjugate gradient method. However, the study does not evaluate the algorithm on large-scale problems and demonstrate scaling only up to 16 MPI processes. Structured nonconvex optimization of large-scale energy systems using PIPS-NLP was performed in [31,12]. The parallel interior point optimization solver for nonlinear programming leverages the dual-block angular structure specific to the problem formulation by applying the Schur complement for efficient parallelization of the linear solves. It was illustrated how different model structures arise in power system domains and how these can be exploited to achieve high computational efficiency. Stochastic optimization problems on high-performance computers have been treated similarly in [30,29]. The numerical experiments suggest that supercomputers can be efficiently used to solve power grid optimization problems with thousands of scenarios under the strict time requirements of power grid operators, particularly due to improved linear algebra on a shared memory level. The proposed additional Schur complement scheme in this work differentiates the proposed formulation from the one suggested by the previous work. The benefits come from a more efficient direct sparse approach for the solution of the underlying sparse KKT system solved at each iteration of the IPM. The proposed method employs one additional Schur complement carefully chosen so that the sparsity of the KKT matrix is maintained. This way the size of the linear system decreases and this allows for memory savings and increased computational performance.

Reduced Space IPM and Adjoint Formulation
The optimization problems in the PDE-constrained problems (e.g., maximization of the oil production [23]) are often solved by reduced space optimization methods [1,3]. The discretized state variables x and the PDEs representing the nonlinear equality constraints (7b) are removed from the optimization problem and are treated explicitly during the evaluation of the objective function value and its gradient. Given the control variables u , the system of the removed equality constraints is used to solve for the state variables x (u ). Thus, the equality constraints are implicitly satisfied. The efficient evaluation of the gradient information is achieved using the adjoint method. The computational cost of evaluating the inequality constraints gradients (7c) can be reduced using the constraint lumping techniques [23]. The second order derivatives are usually not evaluated exactly due to the excessive computational cost, only approximations such as BFGS are used. The concepts of the reduced space IPM are first applied to the OPF problem, the extension to the SCOPF problems is discussed later.
The OPF problem can be interpreted in the context of PDE-constrained optimization problems, where the PF equations (7b) represent the role of the PDEs. The control variables u and the state variables x (u ) correspond to the known and unknown quantities in the PF problem. The known global control variables u contain real power produced at the generators (except the slack generator) and the bus voltage magnitude involved in reactive power balance. The local variables x (u ), representing the state of the power grid, are the remaining quantities such as voltage angles, voltage magnitudes at PQ buses, generator reactive power injections, etc.
The main advantage of the reduced space IPM is that the computational complexity mainly depends on the degrees of freedom, that is, the number of control variables n u = |u |. The reduced space technique is very suitable for solving the optimization problems where the number of control variables n x = |x | is much smaller that the total number of variables, which is the case especially for the SCOPF problem. Furthermore, the reduced space approach is suitable for applications such as time critical optimization, where one may want to terminate an algorithm before optimality has been reached and be assured that the current approximate solution is feasible. In other applications, the objective function might not be defined outside of the feasible region. The main drawback is that in the line search framework, the objective function might be required to be evaluated multiple times with the updated control variables, which requires to resolve the dependent state variables repeatedly. In the OPF and SCOPF problems, such computation is expensive but the evaluation of the gradients will require much more computational resources.
In the reduced space approach, the optimizer controls only u variables but since f (x (u ), u ) and h(x (u ), u ) must be evaluated, the variables x (u ) must be known. The states x (u ) of the problem are obtained as the solution of (7b) for given controls u . Therefore, the network constraints are not forwarded to the optimizer as additional equality constraints, which considerably reduces the size of the optimization problem, not only in terms of the primal state variables but also the associated dual variables for the constraints. The reduced space problem is formulated as

Evaluation of the State Variables
The variables in the full space OPF problems are the complex voltages specified at each bus v = ve j θ and the generator power injections S g = p + j q. We further distinguish three types of buses, namely reference REF, PV and PQ. The known control variables in the reduced space problem are u = (v P V , p P V ), that is, voltage magnitudes at the PV buses and the generator active powers at the PV buses. The unknown state variables . The remaining two variables, θ R E F , v R E F , are fixed and known, therefore are not considered as additional variables.
The unknown state variables x are computed by the standard Newton iterations with a set of the nonlinear mismatch equations g(x ) = 0, defined in (7b), where the control variables u are considered as known, fixed parameters. The equations are ordered as g(x ) = (g 1 (x ), g 2 (x )), where g 1 (x ) = (ℜ(g P V ), ℜ(g P Q ), ℑ(g P Q )) and g 2 (x ) = (ℜ(g R E F ), ℑ(g R E F ), ℑ(g P V )). Considering the aforementioned ordering of the state variables and the mismatch equations, the Newton iterations have the following structure g 11 g 12 g 21 g 22 The block g 12 = ∂ g 1 (x )/∂ x 2 is zero block, since none of the components in x 2 are present in the equations g 1 . The update ∆x 1 can be computed by using the standard power flow problem [35]. If the power grid is connected, the power flow problem g 11 ∆x 1 = −g 1 (x ) is well posed, and the Jacobian g 11 = ∂ g 1 (x )/∂ x 1 is square and non-singular for fixed u . The update ∆x 2 can be computed from

Evaluation of the Gradients
Usage of the gradient-based methods calls for an efficient tool to compute the gradient of the objective function with respect to the design variables, with maximum accuracy and minimum CPU cost. Using the chain rule, the gradient d f d u can be expressed as The term d u x is not known, since x cannot be expressed analytically as a function of u . The solution of (7b) can be obtained using only numerical methods. One method to approximate the gradient is to use n u finite differences over the elements of u in order to approximate Many optimization methods make use also of the second order information, that is, not only the objective gradient but also Hessian is considered. The second order algorithms have a much faster theoretical convergence rate, but this comes with the cost of evaluating the Hessian. An overview of the Hessian evaluation cost is discussed in [28]. Alternatively, one may rely on first order Hessian approximations, e.g. BFGS of DFP [26].

Derivation of the Adjoint Equations
In order to use the gradient based optimization techniques to solve (7), we need to know the first variation of the objective function δ f . The efficient computation is addressed by the adjoint method. It exploits the fact that during the solution process of the state variables (9) the Jacobian matrix ∂ x g is used in Newton's method. The adjoint method uses the transpose of the Jacobian in order to compute the objective function's first variation.
We can introduce the augmented objective function where λ is the vector of Lagrange multipliers. As g(x , u ) is zero by construction, the multipliers may be chosen freely without changing the minimum of the augmented function, therefore, f (x , u ) = (x , u , λ).
Since g(x , u ) is zero everywhere, its gradient is zero as well, and the first variation of f can be expressed as In order to achieve δ = 0 we require ∂ /∂ x = 0 and ∂ /∂ u = 0. To satisfy ∂ /∂ x = 0, we require that the Lagrange multipliers satisfy the following equations: With this choice of the Lagrange multipliers the first term of (14) becomes zero and the gradient of the objective function with respect to the controls that is required by the optimization software is Note that in the process of solving for the state variables (9) the Jacobian of the constraints ∂ g/∂ x is used to determine the Newton step The adjoint equation (15) solves a linear system that differs in form from (17) only by the adjoint operation and the RHS vector. If the power network is connected, the power flow problem g(x , u ) = 0 is square and well posed, and the Jacobian is nonsingular for fixed u .

Detailed Treatment of the Objective Adjoint System
We want to evaluate the objective function gradient (14), where the ordering of the variables and the constraints is defined in section 5. The gradient, considering also the contribution of the augmented function, is Since the objective function depends only on the p variables, there is only a single nonzero entry in the partial derivative w.r.t. the state variables x , and a zero and a dense block in partial derivative w.r.t the control variables

Evaluation of the Multipliers
In order to evaluate the objective function gradient (14), we require that the term in front of δx evaluates to zero by an appropriate choice of the multipliers λ = [λ 1 , λ 2 ]. We require that thus obtaining the adjoint system, The structure of the adjoint linear system (∂ g/∂ x ) , considering ordering of the variables and equations from section 5, and the fact that g 12 = 0, g 22 = −I , is the following We can decouple the solution of λ 1 and λ 2 by solving first for λ 2 , which is trivial. Afterwards, the smaller adjoint system is solved Furthermore, by considering the structure of the RHS vector ∂ f /∂ x 2 , which has only a single nonzero entry at the first position, only the first column of g 21 needs to be evaluated, which corresponds to the gradient of a single constraint ∂ ℜ(g R E F )/∂ x 1 .

Evaluation of the Gradient
Once we have the value of the multipliers, we substitute the multipliers λ in the second term in (14) in order to evaluate the gradient, that is The nonzero pattern of the multipliers is λ 2 = [nn z , 0, 0, . . .], therefore we need to evaluate only the rows of the ∂ g/∂ u corresponding to gradients of constraints [g 1 , ℜ(g R E F )]. Furthermore, ∂ g/∂ p P V is zero, except ∂ ℜ(g P V )/∂ p P V = −I .

Constraint Handling
Constraints that appear as simple bound constraints on the control variables (7e), such as generator power output bounds, can be used directly as inputs to the optimizer. These constraints are removed internally by the optimizer and transformed into standard logarithmic barrier terms. Constraints on state variables (7d), such as the reactive power output, require the solution of the forward PF problem for the evaluation of the constraint and the solution of the adjoint problem for the evaluation of the gradient of the constraint with respect to the control variables. The number of nonlinear inequality constraints (7c) specified in the standard OPF control problem are 2N L , where N L is the number of transmission lines. The inequality constraints consist of two constraints for each line, one for the "from" end and one for the "to" end, imposing limits on apparent power flows expressed in MVA. Computing the gradient for each one of these constraints requires the evaluation of the Lagrange multipliers corresponding to each constraint. Similarly to the equation (12), the augmented constraint i is formulated for each of the inequality constraints, i = 1 . . . 2N L , The first term is required to be zero, therefore the same adjoint system is solved for λ as in (14), but with a different RHS. Since (25) requires the solution of a linear system, 2N L linear systems have to be solved in total to evaluate the gradient of the nonlinear constraints. The bounded state variables are voltage magnitudes at PQ and reference buses, active power at the reference generator, and the reactive power produced at generators. The gradient evaluation of the state variable bounds (7d) is treated similarly as the nonlinear inequality constraints. We define the upper bound constraints as c m a x (x ) = x − x m a x and the first variation of the augmented constraint is Since ∂ c m a x /∂ x is −I , the RHS vectors of the adjoint systems will be columns of the identity matrix. The lower bound constraints c m i n (x ) = x m i n −x are treated in the same way. The total number of adjoint systems solved is 2n x , however, since the lower bound adjoint system varies solely in the sign of the RHS, only n x systems need to be solved.

Detailed Treatment of the Constraints Adjoint System
In order to evaluate the constraint adjoint (25), we follow similar process as in case of the objective function.
That is, we require the term in front of δx to be zero We need to solve the adjoint system, similarly as in the section 6.2. The adjoint system is identical, only the RHS vector vary. The nonlinear constraints h i (x , u ), i = 1 . . . , N L are functions of the θ and v, thus the nonzero entries in the partial derivative w.r.t. the state variables x , and the control variables u are where ∂ h i /∂ x 1 contains nonzero entries and ∂ h i /∂ x 2 = 0. Similarly, ∂ h i /∂ v P V may or may not contain nonzero entries, which is given by the power grid connectivity (e.g. if the PV bus is connected to one, or both ends of the i th transmission line). Finally, ∂ f /∂ p P V = 0. Assuming this nonzero structure, the adjoint system simplifies to The adjoint gradient can be subsequently evaluated substituting the multipliers λ into The nonzero pattern of the multipliers, specifically λ 2 = 0, allows us to evaluate only the rows of the ∂ g/∂ u corresponding to gradients of constraints g 1 and ignore all gradients of constraint g 2 . Furthermore, ∂ g/∂ p P V is zero, except ∂ ℜ(g P V )/∂ p P V = −I .

Constraint Lumping
Explicit evaluation of all inequality constraints (7c) and especially their gradients is computationally intensive for large, realistically sized problems. It suffices for the feasible solution that the maximum of all lines with given rating, l = 1 . . . 2N L , is within the limits, that is, If this maximum is honored, than all of the true constraints are guaranteed to be satisfied. Constraints that are described by nondifferentiable functions can be challenging to incorporate in the optimization frameworks. The approximation of the max function by smooth functions is known in the literature as constraint lumping [27,23]. The max function from (32) is approximated as where α = 0.05h m a x (which needs to be chosen once at the beginning and fixed, not updated in every iteration with new max(h), additionally it needs to be positive) and All the line constraints (7c) in the OPF problem can be thus replaced by a single augmented constraint (33) (assuming that all lines have the same power flow rating, otherwise there will be as many max constraints as there are line categories). The gradient of the constraint, which is required by the optimizer, augmented with the equality constraints (7b), is given by The second term in (35) is required to be zero, therefore, the Lagrange multipliers λ for the constraints are computed from the solution of the adjoint problem Note that only a single adjoint system needs to be solved in order to evaluate the gradient of the lumped line flow constraints, compared to 2N L adjoint systems in (25). It is important to recognize that the approach used for constraint lumping can impact the performance of the optimization procedure. The smaller the coefficient multiplying h m a x in (33), the more accurate the approximation of max becomes. It is important to use a small numerical value, e.g. 0.05, in order to ensure that no overflow occurs in any of the exponential terms in the summation. The approximation of max in (33) is always greater than the maximum of the component line flows, so if this maximum is honored, the true constraint is guaranteed to be satisfied.

Extension to Security Constrained OPF
The formulation of the problem and computation of the adjoint systems in the SCOPF problem need to account for multiple sets of the state variables, corresponding to multiple contingency scenarios. The formulation of the problem (7) The state variables x c (u ) are computed by solving NC+1 PF problems (9). The set of the nonlinear equations g c (x c ) = 0 is slightly different for each scenario, representing the modified grid with the particular contingency c . Note that the computations are independent of each other, therefore can be evaluated in parallel.
The augmented objective function used to formulate the adjoint system is and its first variation is Note the term ∂ f /∂ x c in the equation (39) above. Since the standard SCOPF objective function is defined as the cost of the active power generation, and the active power outputs are the control variables u , the alluded term will be a vector of zeros for all scenarios. Only the active power output of the reference generator is part of the state variables. The reference generator needs to account for different transmission losses in every scenario. If we include the reference generator power output in the cost function only for the nominal case (c = 0), then only a single adjoint system needs to be solved in order to evaluate the gradient of the objective function since the other terms vanish ( ∂ f ∂ x c = 0, therefore, λ c = 0 for c > 0). If the reference generator active power output from every scenario is considered in the objective function, then N C + 1 adjoint systems need to be solved in order to evaluate the gradient of the objective function. For comparison, there was one adjoint system (14) in the OPF problem, which corresponds to considering all generators power output for the nominal case.
The augmented nonlinear line power flow constraint function for the line i in contingency scenario c is defined as (40). A similar process applies to the state variable bounds.
The gradient is expressed as The computational complexity for the SCOPF problem grows as the problem size increases and as more contingency scenarios are considered. Overall, a single adjoint systems needs to be solved for each one of the 2N L line constraints in each contingency scenario c , including the nominal case. The gradient of state variable bounds can be evaluated with the cost of solving n x adjoint systems for each scenario. Overall, for the standard OPF, 1 + 2N L + n x adjoint systems need to be solved; in case of SCOPF there are (N C + 1)(1 + 2N L + n x ) adjoint systems that need to be solved. Introducing the line flow constraints (37b) for individual contingency scenarios may lead to an excessively large number of constraints for realistic networks. This can have an adverse effect for the overall runtime performance of the SCOPF problem. The constraint lumping can be applied, as proposed in the previous section, which reduces the computational burden. The computational time can be further reduced by parallel computation of the gradient information of the constraints associated with each contingency scenario.

Numerical Results
The most expensive component in the full space IPM with the exact Newton is the solution of the large-scale KKT system. While this step is avoided if using the BFGS approximation, the computational complexity of the reduced space quasi-Newton approach is expected to be dominated by evaluation of the state variables and evaluation of the gradients of the objective and constraint functions using the adjoint method, as demonstrated in Figure 1. The performance comparison for pegase9241 benchmark was performed for standard full space method with exact hessian, full space method with the BFGS Hessian approximation and the reduced space method with constraint lumping. The pegase9241 is the largest benchmark available, where the branch power flow limits are specified. The benchmark was run for 3 IPOPT iterations, where the number of line searches was always one, so that the number of function evaluation for all three methods was the same. The precision of the timer is 1ms, so the routines taking less than one millisecond are not shown on the logarithmic plot. Similarly, the OPF Hessian is evaluated only for the exact IPM and the state variables evaluation is performed only in the reduced space method, therefore shown only for the corresponding methods.
Overall, the cost of the three approaches is very similar. The advantage of the reduced space approach is that the expensive components are embarrassingly parallel, thus the solution can greatly benefit from utilizing the parallel computations. However, by dropping the exact Hessian information, the convergence of the IPM might be negatively influenced.

Evaluation of the State Variables
The evaluation of the state variables requires solving the power flow (PF) equations. The PF equations are backbone of the reduced space method. Newton-Raphson (NR) solution algorithm is used, due to its precision and fast convergence. It was observed that the reduced space method is sensitive to the quality of the PF solution, as demonstrated in Figure 2. If the PF tolerance is relaxed, the IPM method needs significantly more iterations until convergence. This possibly rules out computationally less expensive approximations such as Fast Decoupled Load-Flow (FDLF) method, which is based on NR, but greatly reduces its computational cost by means of a decoupling approximation that is valid in most transmission networks (the problem is split into two subproblems with smaller Jacobians which need to be factorized only once). The PF solution in case118 is solved up to the tolerance 1e-14, however for larger networks, such as PEGASE1354, only 5e-12 can be reached. It was also observed that the reduced space method is sensitive to the quality of the limited memory BFGS approximation, as demonstrated in Figure 3. If longer history of the problem is considered, the IPM method needs significantly less iterations until convergence. The reduced space IPM is also sensitive the the barrier parameter update strategy, as demonstrated in Figure 4. The monotone strategy results in significantly more iterations required until convergence.

Constraint Lumping
The computational complexity of the reduced space IPM can be reduced by lumping the constraints representing the transmission lines with the same power flow limit. In the worst case, when each line has a different limit, the complexity would be improved by a factor of two, since there are two constraints for each line (i.e., the two constraints for each line will be replaced by a single constraint, approximating a maximum of the two). The research question is how does the approximation of the maximum influence the convergence of the reduced space IPM.
It is important to note that the state variables (computed as a solution of the PF problem given the control variables) is very sensitive to the perturbation of the control variables. Even for a small perturbation of the controls, the output reactive powers q vary in great degree (while the max perturbation of controls by IPOPT is 2%, the PF solution q change by a factor of 60) resulting in serious q bound violations and related branch power flow heat limits, representing the constraints that are subject to the lumping.
Due to the aforementioned sensitivity of the PF solution to perturbation of the control variables, the output powers q change rapidly. As a consequence, the branch power flows change similarly since they are function of the bus powers. Consequently, the values of the branch power flow violations across individual branches vary significantly, both across individual branches are between the iterations. Three variants of the nonlinear inequality constraints h (x ) and their scaling were investigated, such that h (x ) = SS * −F max L 2 ≤ 0, where S is the complex bus power and the branch power flow limit is denoted as F max L . The original constraints and two scaled versions are the following: In cases 2. and 3. there is one distinguished maximum which is significantly larger from the rest of the values, which makes the lumping approximation more precise. However, the unscaled variant 1. and variant 2. are not appropriate for the lumping, since the large values could easily cause overflow. Therefore, when the large value is encountered the variant 3. is chosen. The disadvantage of variant 3. is that it causes more objective function and inequality constraints evaluations. For this reason, if the maximum value of the scaled constraints (variant 2) is below certain threshold we do not perform additional scaling by maximum as in variant 3. The overall algorithm can be summarized as follows: The benchmarks such as case118, or case2737sop do converge, however PEGASE benchmarks do not. This is true also for the algorithm where no lumping is performed, therefore the convergence problem lies either in the Hessian approximation or the reduced space method itself (maybe insufficient precision of the PF solution).

Full-Space Quasi-Newton Convergence Study
It is well known that quasi-Newton's methods have worse convergence that their exact counterparts. The convergence of the quasi-Newton's full space approach is studied in this section. Figures 5 and 6 illustrate the convergence of both exact and quasi-Newton full-space approaches for two benchmarks. In both cases, the quasi-Newton requires significantly more iterations to reach the required tolerance ε t o l = 10 −2 . However, it is important to note that the quasi-Newton method reaches the optimal value and constraint (primal) feasibility within a similar number of iterations as required by the exact Newton method. The remaining iterations are performed in order to satisfy the tolerance for the dual feasibility (i.e., optimality). We can argue that since we have observed that the objective function is optimal and constraint feasibility is satisfied, the IPM iterations can be stopped early by relaxing the dual infeasibility tolerance. The convergence problems might occur due to the fact that the BFGS cannot approximate the Hessian at late iterates (close to the optimal point) due to the well-known ill-conditioning of the Hessian in the IPM methods as the solution is approached. Exactly for this reason, since we are close to the solution, we can terminate early, not waiting for the dual feasibility reaching the desired tolerance. It is also noted in the IPOPT manual 1 that quasi-Newton methods have trouble bringing down the dual infeasibility. Therefore, special acceptable convergence tolerance is provided. If the usual tolerances are satisfied at an iteration, the algorithm immediately terminates with a success message. On the other hand, if the algorithm encounters many iterations in a row that are considered acceptable, it will terminate before the desired convergence tolerance is met. Such acceptance stopping criterion may be based on the objective function change.

Conclusions
The optimization problems in the PDE-constrained problems (e.g., maximization of the oil production [23]) are often solved by reduced space optimization methods [1,3]. The discretized state variables and the PDEs representing the nonlinear equality constraints are removed from the optimization problem and are treated explicitly during the evaluation of the objective function value and its gradient. Given the control variables, the system of the removed equality constraints is used to solve for the state variables. Thus, the equality constraints are implicitly satisfied. The efficient evaluation of the gradient information is achieved using the adjoint method. The computational cost of evaluating the inequality constraints gradients can be reduced using the constraint lumping techniques [23]. The second order derivatives are usually not evaluated exactly due to the excessive computational cost, only approximations such as BFGS are used.   This approach was applied to the nonconvex OPF problems and the the results demonstrated that the computational cost of evaluating the gradient information is excessive and the constraint lumping introduces non-smooth functions, which leads to convergence difficulties of the IP method used for solution of the OPF problems. Additional investigation needs to be performed, closely analyzing the Hessian approx-imation and smoothing technique used during the constraints evaluation.