Raquel García-Blanco

Pedro Díez
email: pedro.diez@upc.edu

Domenico Borzacchiello

Francisco Chinesta

Algebraic and parametric solvers for the power flow problem: towards real-time and accuracy-guaranteed simulation of electric systems

Keywords: Power Flow problem, Reduced Order Model, Distributed Generation, Optimization, Accuracy control

The power flow model performs the analysis of electric distribution and transmission systems. With this statement at hand, in this work we present a summary of those solvers for the power flow equations, in both algebraic and parametric version. The application of the Alternating Search Direction method to the power flow problem is also detailed. This results in a family of iterative solvers that combined with Proper Generalized Decomposition technique allows to solve the parametric version of the equations. Once the solution is computed using this strategy, analyzing the network state or solving optimization problems, with inclusion of generation in real-time, becomes a straightforward procedure since the parametric solution is available. Complementing this approach, an error strategy is implemented at each step of the iterative solver. Thus, error indicators are used as an stopping criteria controlling the accuracy of the approximation during the construction process. The application of these methods to the model IEEE 57-bus network is taken as a numerical illustration.

Introduction

Power system engineering is a technology field inside the general discipline of energy and electrical engineering. Specifically, power flow analysis is a branch of this area that deals with:

-Transmission and distribution of electrical power -Energy management, storage and generation -Power system planning: operation and expansion -Optimal control and contingency analysis -Real-time monitoring and security risk assessment for reaching stability and reliability -Decision making and voltage regulators-assessment Besides of these aforementioned applications, one of the most significant ones is the design verification. This involves the study of the state and management of the physical networks. The amount of money that should be invested in order to modify, repair and/or expand such networks might be large. According to this, the early simulation of the whole procedure for avoiding unnecessary failures and assuring the viability of the process becomes potentially necessary. Thus, design verification is related to operation and expansion planning.

The design of a network also concerns electricity production, hence terms such as renewable energy and environment emerge. This sort of technology is directly related to another application of the power flow analysis, the optimization of networks under some constrains. Fortunately, renewable energy has a strong impact on the electricity market ensuring the reduction of the greenhouse emissions. This can be done through the introduction of Distributed Generators (DGs) in the case of distribution systems or just generators for transmission systems, for instance, photo-voltaic panels or wind turbines. It is widely believed that the distributed power generation is a technology that could help to enable efficient renewable energy production. DG technology is related to the use of small generating unit installed at strategic points of electric power systems [START_REF] Shrivastava | Optimal Placement Methods of Distributed Generation:A Review[END_REF]. Both the design network and the optimization should be efficient and provide security to guarantee all the desirable benefits. Therefore, these procedures must be analyzed and tested previously. In such way, the Uncertainty Quantification (UQ) measures the error and uncertainties, being another relevant application.

During the last decades, the development of network simulators has been fundamental in the field of power flow problem numerical simulation. However, these improvements have also brought new challenges, for instance, the possibility of making decisions about the state of the grid in real time. This requires solving different configurations of the same problem quickly. The work of this paper focuses on summarizing our contributions to this issue. Our goal is to find the optimal location and sizing of a generator which is set in the network minimizing the system losses in real time. Moreover, we control the quality of the solution in terms of the power losses during its construction.

The layout of the paper is structured as follows: section 2 presents the governing equations distinguishing the algebraic version from the parametric version of the power flow problem while in section 3, a brief review of the deterministic and classical solvers and the Alternating Search Directions method applied to the power flow equation are presented. In this section, the error assessment for the algebraic version of the problem is also illustrated. Equivalently, section 4 analyses probabilistic methods as well as the error assessment associated with the parametric version of the equation. Section 5 details numerical examples where the proposed method was applied. The paper closes with section 6 in which we present our main outcomes regarding the techniques for solving the power flow problem.

Problem Statement

Algebraic problem

The basic formulation of the well-known power flow problem was originally illustrated in [START_REF] Elgerd | Electric Energy Systems Theory[END_REF][START_REF] Stagg | El-Abiad: Computer Methods in Power System-Analysis[END_REF]117]. The main objective of the power flow solution is described by [117] as: to obtain the individual phase voltages at all nodes in the network corresponding to specified system conditions. Consequently, the unknowns of the problem are the voltages and nodal intensities collecting in vectors of n components V and I ∈ C n , where n is the number of degrees of freedom. Note that for a threephase distribution system, n is three times the number of buses. The input data characterizing the power flow problem is the following:

-The topology of the grid, described by the number of lines, the number of buses and their connectivity. -The complex power source vector S ∈ C n , describing the power supplied and/or extracted at each phase of each node. -The admittance matrix Y ∈ C n×n including the material characteristics of the devices conforming the grid. -The vector I 0 ∈ C n accounting for the current originated by the slack node. Introducing a slack node is necessary to guarantee the solvability of the problem. The complex voltage in this node is known, and therefore it is not reevaluated. This is equivalent to reduce the dimension of the admittance matrix by deleting the slack bus row and column, see [START_REF] Dimitrovski | Slack bus treatment in load flow solutions with uncertain nodal powers[END_REF][START_REF] Grainer | [END_REF][START_REF] Kabir | Inclusion of slack bus in newton raphson load flow study[END_REF].

At each bus, the nonlinear relation between the voltage, the current and the complex power is provided by the following equation:

S = V I * , (1)
where I * denotes the complex conjugate of the current vector I, and the symbol denotes the Hadamard product of vectors (component-wise product). Moreover, Kirchhoff's law leads to the following algebraic system of equations:

YV = I + I 0 , (2)
which, using [START_REF]Power systems test case archive[END_REF] results in a nonlinear algebraic system of equations for the unknown V :

YV = S * V * + I 0 = I bus (V ) , (3)
where the symbol denotes the component-wise quotient between vectors.

The admittance matrix and power source in Cartesian form are Y = (Y) + iI(Y) and S = (S) + iI(S) = P + iQ respectively where (•) and I(•) stand for the real and the imaginary part of the matrix or vector and i is the imaginary unit. Thus, the vector of voltages reads V = (V ) + iI(V ). We also adopt the notation V = V Re + iV Im , to shorten some expressions in the following. Moreover, the vector V is also expressed in polar form (module-argument form, α l being the argument of V l ), such that each component reads V l = |V l | (cos(α l ) + i sin(α l )), for l = 1, . . . , n, now the power flow equations (3) read as

                     P l = n k |V l ||V k |[Y Re lk cos(α l -α k )+ + Y Im lk sin(α l -α k )] -P 0 Q l = n k |V l ||V k |[-Y Im lk cos(α l -α k )+ + Y Re lk sin(α l -α k )] -Q 0 , (4)
where the module and argument of the slack node are known, thus the term V * I 0 is also known and from now on it is called S 0 = P 0 +iQ 0 . Furthermore, θ lk = α l -α k is defined as the difference in voltage angle between the l-th and k-th buses, hence the system of equations now is

           P l = n k |V l ||V k |[Y Re lk cos(θ lk ) + Y Im lk sin(θ lk )] -P 0 Q l = n k |V l ||V k |[-Y Im lk cos(θ lk ) + Y Re lk sin θ lk ] -Q 0 . (5)
This is a nonlinear real system of 2n equations and 2n unknowns. For each node in the network, that is l = 1, . . . , n, the active power P l and the reactive power Q l are known, while the |V l | and α l are unknown variables.

Equations ( 3) and ( 5) are equivalent, the choice of one rather than the other depends on the type of nodes in the network and the available data. If there are PQ nodes, where the values of P and Q are known, both equations are used. However, in the case of the PV nodes where just P and |V | are given, equation ( 5) is more suitable.

Parametric problem

The analysis of electric networks under different configurations of loads requires the solution of similar problems a large number of times. Describing these scenarios is easily done by introducing parameters and therefore defining the concept of Parametric Power Flow problem. The novelty of this methodology is the fact that diverse parameters of the power flow problem are now considered as variables rather than input quantities.

Taking into consideration the optimization problem previously described in section 1, typical examples of parameters are the location and nominal power, denoted by q and r, of some generator, and the time t that modulates the power S. In this work, the general form of the Parametric Power Flow problem is described by taking all the variables depending on these three parameters in [START_REF] Acharya | An analytical approach for DG allocation in primary distribution network[END_REF], for instance S(q, r, t). Consequently, its solution also depends on these parameters, namely V (q, r, t). In practice, this brings the problem from a simple nonlinear algebraic equation in C n into a multidimensional setup, further details are shown in section 4.2.

Algebraic solvers

This section summarizes the state-of-art of the most significant algebraic methods for solving the power flow equations. Moreover, the Alternating Search Directions Methods (ASDM) along with a review of the error approach in this field is also introduced.

Y-matrix and Z-matrix methods

Over the last decades, numerous methods have been proposed in order to solve the power flow equation. The first practical technique emerged in 1956 [116]. During the same decade, methods called Y-matrix also appeared [START_REF] Brown | Digital solutions for large power networks[END_REF][START_REF] Glimn | Automatic calculation of load flows. Power apparatus and systems[END_REF]. These methods are a straightforward fixed-point iteration from (3). Thus, an approximated value V [γ] is used to compute the next iteration V [γ+1] such that

YV [γ+1] = [S * V [γ] * + I 0 ] , (6)
by solving, in each iterative step, a linear system of equations with matrix Y. These methods are consistent and were successfully employed in many examples, however, they do not guarantee converge. This difficulty was overcome through the introduction of the Z-matrix methods [14-17, 53, 55]. The main idea of these methods is to invert the system admittance matrix Y, obtaining the impedance matrix Z, using a technique based on Kron's concept of network tearing using the system data. This procedure is faster than the standard matrix inversion and avoids the necessity of complete re-inversion when such minor changes in the network are required. These changes are made directly to the inverted matrix, thus the computation time involved for such modifications is a small fraction of that needed for a complete matrix inversion [START_REF] Brameller | Some improved methods for digital network analysis[END_REF]. Besides, when it comes to networks under fault conditions, the Y-matrix approach requires an iterative solution of the entire network for each fault condition. However, one of the distinct advantages of this method is that (once the matrix is formed) all fault calculations may be obtained with a minimum of arithmetic operations involving only related portions of the matrix [START_REF] Brown | Digital calculation of 3-phase short circuits by matrix method[END_REF].

Many techniques have been proposed to modify the traditional Z-matrix building algorithms. Among those methods, the Gauss implicit Z-matrix method is the most generally used method [START_REF] Ou | A novel Z-matrix algorithm for distribution power flow solution[END_REF]. Furthermore, some novel studies about the convergence analysis of this method including PV nodes and DG have emerged in the current decade [START_REF] Chiang | Convergence/divergence analysis of implicit Zbus power flow for general distribution networks[END_REF][START_REF] He | An improved power flow algorithm for distri-bution networks based on Z-bus algorithm and forward/backward sweep method[END_REF]122,129].

Gauss-Seidel and Newton-Raphson methods

Around the sixties, the notable Gauss-Seidel (GS) [START_REF] Glimn | Automatic calculation of load flows. Power apparatus and systems[END_REF]106] and Newton-Raphson (NR) [START_REF] Ness | Iteration methods for digital load flow studies. Power Apparatus and Systems[END_REF]110] methods for power flow calculations were also proposed. The equations for both iterative methods are:

-Gauss-Seidel At iteration γ, the solution V [γ+1] is obtained solving the below system of equations:

Y L V [γ+1] = S * V [γ] * -Y U V [γ] + I 0 , (7)
being Y U the upper triangular part and Y L the lower triangular part of Y plus its diagonal. In the literature, Gauss-Seidel method is also classified as an Y-matrix method since the calculation of the solution depends on the admittance matrix. -Newton-Raphson Before defining NR equations, since the conjugate function is not a holomorphic function and complex derivation is not formally defined, it is necessary to consider the Cartesian representation of the vectors and matrices involved in the power flow equations. By introducing the real and imaginary parts of voltages, currents, powers and admittances as separate variables, the following quantities are defined as:

Ŷ = Y Re -Y Im Y Im Y Re , V = V Re V Im , Îbus = I Re bus I Im bus = I Re 0 I Im 0 + (P V Re + Q V Im ) (V Re V Re + V Im V Im ) (P V Im -Q V Re ) (V Re V Re + V Im V Im ) , (8)
where Ŷ ∈ R 2n×2n and V , Îbus ∈ R 2n are duplicating dimensions of the complex matrices and vectors. The power flow equations can be written now as:

Ŷ V = Îbus ( V ) . (9)
Newton-Raphson method consists in iteratively updating V with an increment ∆ V , that is

V [γ+1] = V [γ] + ∆ V .
In the following the dependence on γ is eliminated in the superscript to simplify notation. Thus, the resulting algorithm reads

Ĵ∆ V = -Ŷ V + Îbus ( V ) , (10)
where the Jacobian Ĵ is the partial derivative of the right-hand-side of (10) (the residual) with respect to V , that is

Ĵ = Ŷ -Ĵ11 Ĵ12 Ĵ21 Ĵ22 (11)
where Ĵkh for k, h = 1, 2 are diagonal matrices in R n×n (component h of Îbus ( V ) depends only on component h of V , see [START_REF] Atwa | Probabilistic approach for optimal allocation of wind-based distributed generation in distribution systems[END_REF]) such that

Ĵ11 = ∂I Re bus ∂V Re , Ĵ12 = ∂I Re bus ∂V Im , Ĵ21 = ∂I Im bus ∂V Re and Ĵ22 = ∂I Im bus ∂V Im .
Particularly,

Ĵ11 ll = P l (|V l | 2 -V Re l ) + 2Q l V Re l V Im l |V l | 4 , Ĵ12 ll = Q l |V l | 2 -P l V Re l + 2Q l (V Im l ) 2 |V l | 4 , Ĵ21 ll = -Q l |V l | 2 -P l V Im l -2Q l (V Re l ) 2 |V l | 4 , Ĵ22 ll = P l (|V l | 2 -V Im l ) -2Q l V Re l V Im l |V l | 4 , for l = 1, . . . n.
Both GS and NR methods enjoy of low memory usage and competent ratios of convergence, better in the case of NR which has an optimal quadratic rate, although the computational time increases because of the assembling of the Jacobian matrix at every single iteration. For that reason, diverse approaches have appeared over the years. In the case of NR methods, it is worth mentioning different modifications of the original problem introducing decomposition of the Jacobian matrix [START_REF] Cheng | A modified Newton method for radial distribution system power flow analysis[END_REF][START_REF] Garcia | Three-phase power flow calculations using the current injection method[END_REF], reformulating the original equations to accommodate the introduction of generation devices [START_REF] Li | A modified Newton-Raphson power flow method considering wind power[END_REF][START_REF] Sameni | A modified Newton-Raphson method for unbalanced distribution systems[END_REF] or decreasing the computational time thorough the application of third-, fourth-and fifthorder Newton-like methods [START_REF] Derakhshandeh | Application of high-order Newton-like methods to solve power flow equations[END_REF]. Similarly, the GS approach has been improved using block version of its initial equations [START_REF] Maffei | A colored Gauss-Seidel approach for the distributed network flow problem[END_REF]107] or combined with the implicit Z-matrix bus method for unbalanced distribution networks [122]. The initial references of these three type of methods (Z-matrix bus, GS and NR) are shown in the reviews [START_REF] Laughton | Numerical techniques in solution of power-system load-flow problems[END_REF]100] while an extensively recent study of them can also be found in [START_REF] Gómez-Expósito | Electric energy systems: analysis and operation[END_REF].

Additionally to the improved methods mentioned above, since NR methods emerged, a great effort has been made to overcome the problem of updating the Jacobian when the size of the test systems was considerable large. Consequently, a variety of formulations have been developed. These include:

-Newton-Krylov methods consisting in solving the Jacobian equation partially or combined with a Krylov subspace method [START_REF] Idema | Fast newton load flow[END_REF][START_REF] Idema | Scalable Newton-Krylov solver for very large power flow problems[END_REF][START_REF] Idema | Towards faster solution of large power flow problems[END_REF]123] -Jacobian-free Newton-Krylov methods (JFNK) where a Krylov subspace is built up for correcting the Jacobian in NR strategy [67] -Jacobian-free methods that analyze approaches as partial Jacobian update variants and inexact solutions [START_REF] Chen | A Jacobian-free Newton-GMRES (m) method with adaptive preconditioner and its application for power flow calculations[END_REF][START_REF] De Leon | Iterative solvers in the Newton power flow problem: preconditioners, inexact solutions, and partial jacobian updates[END_REF] Generally, despite the fact that these strategies does not include the whole Jacobian, the quadratic convergence is still granted.

Apart from these approaches, the most popular is Fast Decoupled Load Flow Method (FDLF) [101]. It mainly consists in approximating the Jacobian by using factorization, preconditioners or information obtained from the first iteration in order to quickly solve the Jacobian system. In such a way, the matrices are kept constant hence NR method is reduced to a sequence of decoupled linear problems for the voltage magnitude and phase angle.

The application of the NR method for the equation (5) results in a nonlinear real system of 2n equations with 2n unknowns, the vectors |V l | and α l for l = 1, . . . , n. The equations read as:

J ∆α ∆|V | = - P res Q res , (12)
where P res and Q res are the residuals of the equation and the Jacobian is

J = J 11 J 12 J 21 J 22 = ∂Pres ∂α ∂Pres ∂|V | ∂Qres ∂α ∂Qres ∂|V | . (13)
Specifically,

[J 11 ] lk =        |V l ||V k |[Y Re lk sin(θ lk ) -Y Im lk cos(θ lk )], k = l - m =l |V l ||V m |[Y Re lm sin(θ lm ) -Y Im lm cos(θ lm )], k = l [J 12 ] lk =        |V l |[Y Re lk cos(θ lk ) + Y Im lk sin(θ lk )], k = l m =l |V m |[Y Re lm cos(θ lm ) + Y Im lm sin(θ lm )] + 2Y Re ll |V l |, k = l [J 21 ] lk =        -|V l ||V k |[Y Re lk cos(θ lk ) + Y Im lk sin(θ lk )], k = l - m =l |V l ||V m |[Y Re lm cos(θ lm ) + Y Im lm sin(θ lm )], k = l [J 22 ] lk =        |V l |[Y Re lk sin(θ lk ) -Y Im lk cos(θ lk )], k = l m =l |V m |[Y Re lm sin(θ lm ) -Y Im lm cos(θ lm )]+ -2Y Im ll |V l |, k = l for l, k, m = 1, . . . n.
Taking into account that for k = l,

[J 11 ] lk = |V k | [J 22 ] lk and [J 21 ] lk = -|V k | [J 12 ] lk ,
the algorithm is rewritten as

J ∆α ∆|V | |V | = - P res Q res , (14)
where

J = J11 J12 J21 J22 , explicitly, (15)
J11 lk = [J 11 ] lk , k = l -Q l + Y Im ll |V l | 2 , k = l J12 lk = -[J 21 ] lk , k = l P l + Y Re ll |V l | 2 , k = l J21 lk = [J 21 ] lk , k = l P l -Y Re ll |V l | 2 , k = l J22 lk = [J 11 ] lk , k = l Q l -Y Im ll |V l | 2 , k = l for l, k = 1, . . . n.
As mentioned before, FDLF method is a variation of Newton-Raphson method. It is achieved by only inverting the Jacobian matrix once it is simplified assuming the below statements:

-It was observed that real power P was barely influenced by changes in voltage magnitude V , thus, all the derivative are considered to be zero. Similarly, Q was relatively insensitive to changes in α. This means that J12

lk = J21 lk = 0.
-The difference between angles θ lk = α lα k is usually small, cos(θ lk ) is taken by 1 and sin(θ lk ) as 0, for l, k = 1, . . . , n. -The magnitude of some voltages is also assumed to be 1.

Applying these assumptions to equation ( 14) and dividing equations by |V k | in both sides, the system to be solved is:

U 0 0 U ∆α ∆|V | = - Pres |V | Qres |V | , (16)
where U = -Y Im and U is built taking the elements of -Y Im that correspond to the PV nodes.

Later, FDLF has been developed for unbalanced radial distribution system [131], three phase distribution networks [START_REF] Lin | Three-phase distribution network fast-decoupled power flow solutions[END_REF] and for transmission system using an optimal multiplier [START_REF] Bijwe | Robust three phase fast decoupled power flow[END_REF]. From the mathematical point of view, some authors have addressed its theoretical background [START_REF] Monticelli | Fast decoupled load flow: hypothesis, derivations, and testing[END_REF]120].

Holomorphic Embedding Load Flow methods

Besides this deficiency in terms of the Jacobian assembly, the traditional power flow methods may suffer from the fact that there is no guarantee that they converge to the physical or high voltage solution. Some iterative solvers might converge to spurious non operative solutions or simply fail to converge in a number of cases. The reason behind such behavior could be either the dependency between the initial estimate and the final approximation [START_REF] Iwamoto | A load flow calculation method for ill-conditioned power systems. Power Apparatus and Systems[END_REF][START_REF] Schaffer | A nondiverging polar-form Newton-based power flow[END_REF][START_REF] Stott | Effective starting process for Newton-Raphson load flows[END_REF] or the system operability making the algorithm not able to find the operative solution. This may happen when the value of network parameters move outside of the standard operating range due to contingencies [114]. In the case of NR methods, the authors of [108,109] showed that the nature of the power flow solution is fractal.

Overcoming both adversities was a challenge which motivated numerous authors. On one hand, methods based on truncated Taylor expansions in a polar or Cartesian coordinate form were proposed [START_REF] Sauer | Explicit load flow series and functions[END_REF][START_REF] De Souza | Non-iterative load-flow method as a tool for voltage stability studies[END_REF]121]. A suitable one is the second order load flow technique which requires less iterations and have better convergence characteristics than conventional NR technique [START_REF] Sachdev | A second order load flow technique[END_REF]. On the other hand, it was recently introduced the Holomorphic Embedding Load Flow Method (HELM) [111,112]. HELM is based on an analytical continuation, a technique that extends the domain of analytic functions, from complex analysis relying on Padé approximants. The method extends the voltage variables into analytic functions in the complex plane providing a non-iterative procedure for constructing the complex power series of voltages.

If a simple two buses system is considered, the Zmatrix method applied to the scalar version of equation (3) reads as:

V = V 0 + Z[S * V * ], (17)
where V 0 = ZI 0 . Rewriting equation ( 17) using the notation U = V /V 0 , the following equation is obtained,

U = 1 + σ U * , (18)
where σ = ZS * |V0| 2 . This above equation ( 18) can be seen as a continued fraction approximation of the solution,

U = 1 + σ 1 + σ * 1 + σ 1 + ... . ( 19
)
As we previously mentioned, the Holomorphic Embedding method is based on analytical continuation and a continued fraction is defined. In this particular case, [START_REF] Chaturantabut | Application of POD and DEIM on dimension reduction of nonlinear miscible viscous fingering in porous media[END_REF] is also seen as the same continued fraction resulting from the application of the Holomorphic Embedding method to the same system, see [111]. This continued fraction suggests the use of Padé approximants and its convergents corresponds to the application of the fixed point equation [START_REF] Brown | Digital calculation of 3-phase short circuits by matrix method[END_REF]. Therefore, the iterative solutions found with the Z-matrix method, coincide with the ones found with the HELM as the number of coefficients of the Padé approximant is increased. In a general case, Holomorphic Embedding changes σ by sσ in equation ( 18) and defines a system of two equations: [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF] with F (s) = F * (s * ). In this way, the functions F (s) and F (s) are holomorphic. Note that F (s = 1) recovers the solution U of equation [START_REF] Brown | Digital solutions for large power networks[END_REF]. The procedure is to consider the power series expansion of F (s) about s = 0 since F (s) and F (s) are holomorphic. The embedded equations [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF] allow to seek the coefficients of the power series as the solution to a succession of linear systems. Particularly, the derivatives of the function F (s) evaluated in s = 0 are:

F (s) = 1 + sσ F (s) F (s) = 1 + sσ * F (s)
F (0) = 1 F (1) (0) = σ F (2) (0) = -2σσ * F (3) (0) = 6(σ 2 σ * + σ(σ * ) 2 ) F (4) (0) = -72(σ 2 (σ * ) 2 ) -24(σ(σ * ) 3 + σ 3 σ * )
F (5) (0) = -60(σ 2 σ * ) + 600(σ 3 (σ * ) 2 ) + 720(σ 2 (σ * ) 3 )+ + 120(σ 4 σ * ) . . . [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF] Using these derivatives, the Padé approximation P (s) is computed and the solution U = F (s = 1) is approximated by P (s = 1), that is to say, using Padé approximants, the solution at s = 1 can be constructed.

The Padé approximants are a particular type of rational approximation for power series. They have been extensively used since their convergence has been known to be much better than the one of power series. For instance, these approximants are usually superior to Taylor series when the functions to be approximated are complex with singularities (poles), because the use of rational functions allows them to be well-represented. In the case of the power flow equation, Stahl's results reveal that Padé approximants are suitable for analytic continuation. In fact, these results confer the method very strong additional guarantees: if the approximants converge at s = 1, the result is guaranteed to be the analytic continuation of the high voltage branch at s = 1; conversely, if the Padé approximants do not converge at s = 1 then it is guaranteed that there is no solution (that is, the system is beyond voltage collapse). For more details, see [111].

After this initial proposal, an extension from alternating current to direct current-based systems has been presented [113]. Other authors have also explored this approach [START_REF] Rao | The holomorphic embedding method applied to the power-flow problem[END_REF]103]. The main advantage of this sort of strategies is its reliability finding a stable solution for any set of power flow equations. If the starting solution is an operative one, there is guarantee that the algorithm converges fast to a solution which is in the branch of the operative solutions.

Alternating Search Directions Method for the Power Flow equation

Having described the power flow equation as a combination of a nonlinear local constraints (1) and a linear global problem [START_REF] Abido | Optimal power flow using particle swarm optimization[END_REF], the first idea is to consider the system formed by both equations, instead of the primitive formulation [START_REF] Acharya | An analytical approach for DG allocation in primary distribution network[END_REF]. Following this path of dividing the original system into two different equations the Method of Alternating Search Directions emerged. More detailed theoretical background is provided in [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics: New Approaches and Non-Incremental Methods of Calculation[END_REF] where the same method is applied to nonlinear structural mechanics problems.

Mathematical background

The application of this strategy in the power flow field generates an iterative solver with two steps, one per equation. For this, additional linear relations between voltages and currents are needed, the so-called search directions, named as matrices α and β.

At iteration γ, for a given matrix α ∈ C n×n and initial pair (V, I) [γ] , an intermediate solution (denoted by superscript γ + 1 2 ) is found from the linear system

I [γ+ 1 2 ] -I [γ] = α(V [γ+ 1 2 ] -V [γ] ) Y V [γ+ 1 2 ] = I 0 + I [γ+ 1 2 ] . (22)
Similarly, in a second step, for a given diagonal matrix β ∈ C n×n , the solution is updated by solving the system

I [γ+1] -I [γ+ 1 2 ] = β(V [γ+1] -V [γ+ 1 2 ] ) V [γ+1] * I [γ+1] = S * . ( 23
)
There are several advantages associated with this approach:

-The non-linearity and the non-locality can be tackled separately, since subproblem (2) is global but linear, while subproblem (1) is nonlinear but local, meaning that each nodal equation can be solved individually.

-If α is constant, the matrix factorization needed to solve system [START_REF] Chen | State estimation and power flow analysis of power systems[END_REF] is only performed once. -Equations ( 23) can be solved analytically.

-Since the pairs (V, I)

[γ+ 1 2 ] and (V, I)

[γ+1] fulfill equations ( 2) and (1) respectively, the algorithm is numerically consistent.

Further details about the convergence of the method, a geometrical interpretation and also an appendix describing the procedure for seeking the high voltage solution are shown in [START_REF] Borzacchiello | Unified formulation of a family of iterative solvers for power systems analysis[END_REF].

Choices of the search directions

Without losing accuracy, currents can be eliminated from the equations ( 22) and [START_REF] Chen | Probabilistic load flow: A review[END_REF], and the iterative algorithm for V can be formulated as follows:

Method

Choice for α Choice for β

Gauss-Seidel Y U ∞ Newton-Raphson ∂I bus ∂V ∞ Z-matrix bus 0 ∞ -Global
Step. Starting from the iterate V [γ] , the intermediate solution V [γ+ 1 2 ] is found by solving the linear system

(Y -α) V [γ+ 1 2 ] = S * V [γ] * -αV [γ] + I 0 . (24)
-Local Step. The new iteration V [γ+1] is then obtained from the solution of the system

βV [γ+1] + (Y -β) V [γ+ 1 2 ] -I 0 -S * V [γ+1] * = 0 . (25)
Depending on the choice of the matrices α and β it is possible to recover some classical methods as is shown in table 1. In the light of the results presented in [START_REF] Borzacchiello | Unified formulation of a family of iterative solvers for power systems analysis[END_REF], and based on the choices

α = diag S * |V b | 2 β → ∞ , (26)
where V b is the voltage base, the proposed approach is optimal since well designed grids are normally operating not far from this point. Besides this, another argumentative reasons for this affirmation are given in [START_REF] Borzacchiello | Unified formulation of a family of iterative solvers for power systems analysis[END_REF].

Algebraic version of the error assessment

Despite that fact the power flow equations have been studied in detail, the errors during the simulations have not received the corresponding attention. Particularly, the error has been addressed from another points of view as identifying errors associated with power controller parameters [130] or taking into account state estimation method for measurement error and model accuracy [START_REF] Amini | DC power flow estimation utilizing bayesian-based lmmse estimator[END_REF][START_REF] Chen | State estimation and power flow analysis of power systems[END_REF][START_REF] Rouhani | Parzen window density estimator-based probabilistic power flow with correlated uncertainties[END_REF].

When it comes to an error in the application of Reduce Order methods (more information is given in section 4), [START_REF] Rathinam | An iterative method for simulation of large scale[END_REF] provides an error analysis of the computed solution of a reduced model obtained from Proper Orthogonal Decomposition (POD), illustrating the method using a power grid example modeled by nonlinear swing equations. In addition, [START_REF] Galbally | Non-linear model reduction for uncertainty quantification in large-scale inverse problems[END_REF] provides an analysis of the errors involved in solving a nonlinear initial value problem using a POD reduced order model. An error bound on the Discrete Empirical Interpolation Method (DEIM) approximation is provided in [START_REF] Chaturantabut | Discrete empirical interpolation for nonlinear model reduction[END_REF][START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF], while [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF] gives an error analyses for the empirical interpolation procedure and [119] presents a-posteriori error estimation for POD-DEIM reduced nonlinear systems. On the other side, the error estimation in Proper Generalized Decomposition (PGD) is still an open question. In this regard, some strategies have been proposed in [START_REF] Ammar | Coupling finite elements and proper generalized decompositions[END_REF][START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF] .

This section focuses on the development of the error equations based on a Quantity of Interest (QoI), particularly the system losses. The classic strategy that has been applied in this work, has been also applied to different problems in the field of error estimation for Reduced Order Models, see [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Florentin | Adaptive reduced basis strategy based on goal oriented error assessment for stochastic problems[END_REF][START_REF] Mozolevski | Goal-oriented error estimation based on equilibrated-flux reconstruction for finite element approximations of elliptic problems[END_REF].

Having defined the error as the subtraction between the two n components vectors: V the actual solution of the problem and V a an approximation, such as

E = V -V a , (27)
the residual equation of ( 3) associated with the approximation is:

R(V a ) = S * -V * a (YV a -I 0 ). (28)
Assuming that R(V ) = 0, the derivation of the error equation is straightforward. The first step is to linearize R(•) by neglecting the quadratic term. After this, the expanded expression reads as

R(V ) = R(V a ) -V * a YE- E * (YV a -I 0 ) = R(V a ) -AE -BE * , (29)
where A = Diag(V * a )Y and B = Diag(YV a -I 0 ) are matrices in C n×n . The operator Diag(•) is introduced to compact the notation such that it produces a square matrix with the elements of a vector on the diagonal.

Clearly, equation ( 29) is still nonlinear because it involves the conjugate operator. In order to linearize it, vectors and matrices are separated in their Cartesian representation. Hence, equation ( 29) is rewritten as a linear system of 2n real equations and unknowns, namely

C Ê = R(V a ), (30)
where the matrix C ∈ R 2n×2n and the vectors R(V a ) and Ê in R 2n . The objective of the optimization problem in this work is to minimize the power losses, thus, they are assumed as the quantity of interest. In general, the positive number representing the losses associated with a vector V is:

l(V ) = (V * T Y L V ) Re (31)
where Y L is the admittance matrix corresponding to the grid, accounting for all the lines and buses but not including the terms associated with the generators. By following the same procedure, the form l(•) is also nonlinear and has to be linearized in order to define a goal-oriented error assessment strategy. The linearized version of the losses after neglecting the second degree term reads,

l(V ) = l(V a ) + (V * a T Y L E) Re + (E * T Y L V a ) Re = l(V a ) + (f T E) Re + (g T E * ) Re , (32)
where f = Y T L V * a and g = Y L V a are vectors in C n . Now, using the Cartesian representation in equation [START_REF] De Leon | Iterative solvers in the Newton power flow problem: preconditioners, inexact solutions, and partial jacobian updates[END_REF], the linear approximation for the error in the QoI is defined as:

E QoI = l(V ) -l(V a ) ≈ λT Ê ( 33
)
where

λ = f Re + g Re -f Im + g Im ∈ R 2n .
As a standard strategy in the error assessment procedure, the dual or adjoint problem is introduced in order to obtain a representation of the error in the quantity of interest:

C T ρ = λ. ( 34
)
The solution of this problem ρ is a real vector of dimension 2n. Assuming that the approximation (33) holds, using ρ and (30), the error in the QoI is readily represented as:

E QoI = λT Ê = ρT R(V a ). ( 35
)
The authors of [START_REF] García-Blanco | Monitoring a PGD solver for parametric power flow problems with goal-oriented error assessment (article in review)[END_REF] detail the difficulties in the linearization of equation ( 29) and also explains the possibility of calculating the matrix C and the vectors R(V a ), λ and ρ just once during the whole iterative process.

Parametric solvers

So far, the present paper focused algebraic solvers for power flow equations. In this section, the attention is drawn to solvers where a parametric representation of the problem is involved. Two particular cases are the Probabilistic Load Flow (PLF) and the Optimal Power Flow (OPF) where a power flow solver is called as many times as particular system configurations need to be evaluated.

The concept of Probabilistic Load Flow was first proposed in the seventies taking into consideration uncertainty of the nodes data [START_REF] Borkowska | Probabilistic load flow[END_REF]. Another historical reference, where the definition of Stochastic Load Flow appeared for the first time, is [START_REF] Dopazo | Stochastic load flows[END_REF]. Since then, several scientific contributions have been published in this regard, in review [START_REF] Li | Comparison between two probabilistic load flow methods for reliability assessment[END_REF], it is claimed that probabilistic load flow methods can be divided into three categories:

-Simulation methods: the main example of a simulation method is the so-called Monte Carlo method [START_REF] Fang | A stochastic power flow method based on polynomial normal transformation and quasi monte carlo simulation[END_REF], which simulates power flow calculations based on deterministic samples. It is well-known as a flexible and robust method, nevertheless is highly timeconsuming because of the need of repeating calculations.

-Analytical methods: based on convolution techniques [START_REF] Allan | Probabilistic techniques in a.c. load-flow analysis[END_REF] or cumulant method [127] are claimed to be more effective computationally. -Approximate methods: the most common are the method of moments and the point estimate method [102].

Apart from this classification, another remarkable methods for solving the probabilistic load problems using techniques as combinatorics [START_REF] Rau | Solution of probabilistic load flow equations using combinatorics[END_REF] are: Hybrid Latin Hypercube Sampling and Cholesky Decomposition [125], and polynomial normal transformation and Quasi Monte Carlo Simulation [START_REF] Fang | A stochastic power flow method based on polynomial normal transformation and quasi monte carlo simulation[END_REF]. Over the past decades, PLF has been applied to different problems as branch outages, photo-voltaic and wind power through distributed generators, wind farm power generation, energy storage, reliability, distribution system planning-design-analysis, although without doubt the most significant is the optimization.

The Optimal Power Flow solution was presented in the sixties [START_REF] Dommel | Optimal power flow solutions[END_REF]. The idea of the classical OPF is a power flow problem in which certain controllable variables are to be adjusted to minimize an objective function such as the cost of active power generation or losses, see [104]. In fact, [START_REF] Georgilakis | Optimal distributed generation placement in power distribution networks: Models, methods, and future research[END_REF] defines the optimal distributed generation placement problem (ODGP) claiming that it provides the best locations and sizes of DGs to optimize electrical distribution networks. When ODGP is solved, the objective function can be single or multi-objective. The main single-objective functions are: minimization of energy losses, minimization of system average interruption duration index (SAIDI), minimization of cost, etc. On the other hand, in [START_REF] Shareef | A review on models and methods for optimal placement of distributed generation in power distribution system[END_REF], ODGP multi-objective formulations are classified as multi-objective function with weights, goal multi-objective index and multi-objective formulation considering more than one often contrasting objectives. When it comes to minimize the annual losses (as we aim in this work), new methodologies have been proposed based on the optimal allocation of DG units in the distribution system with the same objective [START_REF] Atwa | Probabilistic approach for optimal allocation of wind-based distributed generation in distribution systems[END_REF][START_REF] Griffin | Placement of dispersed generation systems for reduced losses[END_REF][START_REF] Martinez | Optimum placement of distributed generation in three-phase distribution systems with time varying load using a Monte Carlo approach[END_REF].

Improved versions of the OPF emerged later on [124] and with this also the number of applications increase notably. Some of them are based on the same ideas as the Probabilistic Load Flow applications, other are economic and pollution dispatch or maximum interchange, however one stands out the optimization problem in presence of distributed generation, mainly wind farms or turbines. Diverse strategies has been addressed in order to solve this particular application as it can be seen in [START_REF] Shrivastava | Optimal Placement Methods of Distributed Generation:A Review[END_REF] by using the following techniques:

-Analytical: zero point analysis focusing on the point of the feeder where the power flow is zero or the 2 3 ruled used for capacitor placement in radial distribution system [118] are examples of this category. -Exact formulas: references [START_REF] Acharya | An analytical approach for DG allocation in primary distribution network[END_REF]88] present exact methods such as the exact loss formula or the gradient method. -Evolutionary: Monte Carlo method [START_REF] El-Khattam | Stochastic power flow analysis of electrical distributed generation systems[END_REF], Hereford Ranch Algorithm (HRA) and Genetic Algorithm (GA) based on genetic concepts [START_REF] Kim | Dispersed Generation Planning Using Improved Hereford Ranch Algorithm[END_REF][START_REF] Mithulananthan | Distributed Generator in Power Distribution Placement System Using Genetic Algorithm to Reduce Losses[END_REF], Simulated Annealing (SA) a local search algorithm [START_REF] Gandomkar | A combination of genetic algorithm and simulated annealing for optimal dg allocation in distribution networks[END_REF], Fuzzy System algorithm built using the fuzzy set theory [START_REF] Kim | Dispersed generator placement using fuzzyga in distribution systems[END_REF], Ant Colony optimization specially designed to deal with large search spaces since it dynamically creates the search routes such as real ants do [START_REF] Favuzza | Optimal electrical distribution systems reinforcement planning using gas micro turbines by dynamic ant colony search algorithm[END_REF], Tabu Search that explores the whole solution space randomly based on the local search [80] and Particle Swarm optimization inspired by social behavior of bird flocking among others [START_REF] Abido | Optimal power flow using particle swarm optimization[END_REF].

The optimal placement methods of distributed generators can be solved using the probabilistic approaches mentioned above or the deterministic ones reviewed in section 3. Further information about methods and techniques proposed for solving OPF and PLF are shown in [START_REF] Chen | Probabilistic load flow: A review[END_REF][START_REF] Frank | Optimal power flow: a bibliographic survey I[END_REF][START_REF] Huneault | A survey of the optimal power flow literature[END_REF]].

Standard Reduced Order Model methods

A parametric version of the power flow equation ( 3) is presented in section 2.2. Having this vision of the problem implies that the solution of optimization or uncertainty quantification among others problems in real-time is reachable. However, in practice, the solution of high dimensional problems may become complicated due to the exponential increase of the degrees of freedom. These type of methods are potentially subject to the curse of dimensionality, that is, to a dramatic increment of the computational cost with the number of dimensions. In D dimensions if each parameter assumes d possible states, the extensive exploration of the parametric space is associated to a volume of information that scales with d D . In this context, Reduced Order Models (ROM) are especially indicated to remedy this deficiency.

A historical review of ROM is given by [START_REF] Pinnau | Model Order Reduction: Theory, Research Aspects and Applications[END_REF]. Originally, they were developed in the area of systems and control theory that studies properties of dynamical systems. The fundamental methods in ROM area were published in the eighties and nineties of the last century.

Proper Orthogonal Decomposition method, also known as Karhunen-Loève decomposition or principal component analysis, was proposed by [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. I-III[END_REF]. The basis theory of this strategy is shown in [START_REF] Pinnau | Model Order Reduction: Theory, Research Aspects and Applications[END_REF][START_REF] Rios | A multidimensional residual functional for obtaining the Proper Orthogonal Decomposition coefficients in model reduction[END_REF]. This method essentially supplies an orthonormal basis for representing the given data in a optimal sense, that is to say, given V1 , . . . , Vnsn ∈ R 2n vectors, we can approximate the so-

lution voltage V by V ≈ V P OD = nsn m=1 ω m Vm = Mω (36)
where M = [ V1 , . . . , Vnsn ] is a matrix and the vector ω ∈ R nsn is the new unknown (instead of V ). In practical application, POD methods make essential use of empirical data taken from numerical simulation. Hence, considering a given set of pre-computed voltage solutions V1 , . . . , Vnsn , called snapshots, the matrix of snapshots is defined as V sn = [ V1 , . . . , Vnsn ] ∈ R n×nsn and used instead of the matrix M in equation [START_REF] Dopazo | Stochastic load flows[END_REF], where n sn is the number of snapshots considered. The choice of the data set plays a crucial role and relies either on intuition or simulations. It is affirmed that the incorporation of empiric data of the original model is one of the advantages of the POD method. POD falls into the category of projection methods where the system is projected onto a subspace of the original phase space. Later on, singular value decomposition (SVD) of the snapshot matrix V ns is carried out in order to obtain an approximation of a set of orthogonal basis functions spanning the solution space, that is to say,

V sn = PΣQ T ( 37
)
where P = [u 1 , u 2 , . . . , u n ] is a matrix with orthonormal columns, Q = [q 1 , q 2 , . . . , q n ] is a orthogonal matrix containing the singular vectors, and Σ is a diagonal matrix with positive or zero elements λ l , l = 1, . . . , n sn , called the singular values. In general, rank(V sn ) min(n, n sn ), so each Vh , h = 1, . . . , n sn can be written as:

Vh = rank(V sn ) m=1 ω m Vm ≈ V P OD = n RB m=1 ω m Vm = Pω (38)
where n RB << rank(V sn ) selected such that accumulative variance

n RB l=1 λ l / nsn l=1
λ l is close to one and, P = [ V1 , . . . , Vn RB ] is a n × n RB matrix and the vector ω ∈ R n RB is unknown.

The connection between POD and SVD lies in the fact that the approximating POD basis should contain as much information as possible. Mathematically, the problem of approximating the snapshot vectors by a single vector is written as a constrained optimization problem. Using the Lagrangian formalism, we derive that a necessary condition for this problem is given by the eigenvalue problem. The singular value analysis yields that V1 solves this eigenvalue problem and the functional value is indeed λ 1 . Now, we iterate this procedure and by construction it is clear that for every n RB ≤ n the approximation of the columns of V sn by the first n RB singular vectors is optimal in the leastsquares sense. Altogether, this leads the way to the practical determination of a POD basis of rank n RB .

This method have been developed in different areas, some of them are: image processing, data compression, signal analysis, modeling and control of chemical reaction systems, turbulence models, control of fluids, electrical power grids, pattern recognition, wind engineering, etc. Despise POD methods are useful in many cases, when non-linear systems are involved, difficulties might appear since the cost of evaluating the smaller system resulting after apply the method still depends on the number of variables of the full model. For this reason, other methods have been developed in recent decades, for instance Discrete Empirical Interpolation Method, meant to be an improvement of the POD approximation because (based on a projection combined with interpolation) it achieves a reduction of the nonlinear term with a complexity proportional to the number of reduced variables [START_REF] Chaturantabut | Application of POD and DEIM on dimension reduction of nonlinear miscible viscous fingering in porous media[END_REF].

In general, the concept of ROM has been known for long time in the field of power system engineering. Indeed, grid equivalencing techniques like Ward reduction [115] or POD [START_REF] Parrilo | Model reduction for analysis of cascading failures in power systems[END_REF] are commonly used to reduce the computational cost of power flow analysis of large systems. The combination of both POD and DEIM methods has been applied for model order reduction for semiconductors in electrical networks using DEIM to treat the reduction of nonlinear components [START_REF] Hinze | Discrete empirical interpolation in POD model order reduction of drift-diffusion equations in electrical networks[END_REF]. Also electrical, thermal and micro-electromechanical systems have been also studied [START_REF] Hochman | A stabilized discrete empirical interpolation method for model reduction of electrical, thermal, and microelectromechanical systems[END_REF]. More recently, works dealing with either OPF or PLF using order reduction techniques rely on Sparse Grid approaches [START_REF] Lin | Uncertainty quantification in dynamic simulations of large-scale power system models using the highorder probabilistic collocation method on sparse grids[END_REF]105,126] or Sparse Tensor Recovery [128]. Both techniques can be classified as collocation approaches, since the solution is reconstructed in the high dimensional space from the values it assumes in a set of particular and well-chosen points called the collocation points.

In this work, ROM techniques are not intended to reduce the degrees of freedom of the physical system but the computational complexity associated to the resolution of high-dimensional parametric equations. For that reason, Proper Generalized Decomposition technique is suitable for the power flow problem. In [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF], it is claimed that the technique is based on DEIM, and thus the nonlinear term is interpolated using the reduced basis instead of being fully evaluated. It is remarkable that although PDG is based on these methods, it is a priori model reduction method because it does not depend on previously computed snapshots. PGD discovers the true dimensionality of the model as a part of the solution of the parametrized equations and does not need train simulations. A careful treatment of the nonlinearity is required when using PGD. Diversified strategies exist, depending on the problem at hand, and can be found in the specialized literature [START_REF] Chinesta | On the use of Proper Generalized Decompositions for multidimensional models[END_REF][START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF]. In this work, we combined the non-linear algebraic solver illustrated in section 3.4 with the PGD strategy.

Proper Generalized Decomposition for the Power Flow problem

The goal of this section is to apply the PGD technique to the power flow equations. The output of PGD is a full parametric solution in a compact separated-variables format. As has been mentioned before, the objective is to calculate the optimal location and sizing of a generator set in the network that minimizes the system losses. Thus, the parameters are the location and the nominal power of the generator, denoted by q and r respectively, and also the time t accounting for the hours in a year. The solution of equation ( 3) is sought using the iterative algorithm explained in section 3.4 combined with a PGD approximation.

It is assumed that the input data S(q, r, t) is represented as a separated variables,

S(q, r, t) = H h=1 α h S S h Qh (q) Ȓh (r) T h (t) , ( 39
)
where H is the number of terms in the S expansion, and for h = 1, . . . , H, α h S are positive scalars, S h ∈ C n are the unit vector modes of powers, and Qh (q), Ȓh (r) and T h (t) are the unit parametric modes. We also assume that the PGD approximation V a has a separated form, that means, V a is a sum of M terms, each of them being the product of functions depending on only one of the parameters, namely

V (q, r, t) ≈ V a (q, r, t) = M m=1 α m V V m Q m (q)R m (r)T m (t) , (40)
where, for m = 1, . . . , M , α m V are positive scalars, V m ∈ C n are the unit vector modes of voltages, and Q m (q), R m (r) and T m (t) are the unit parametric modes. The modes are normalized and the positive scalar α m V collects the amplitude of each term.

In practice, the parametric dimensions are discretized in a Finite Element fashion. Let n q , n r and n t denote the number of degrees of freedom discretizing the three parametric dimensions. Thus, function Q m (q) is identified with vector Q m ∈ C nq , similarly vectors R m ∈ C nr and T m ∈ C nt represent functions R m (r) and T m (t). Thus, the multivariate function V a (q, r, t) (from I q ×I r ×I t to C nr ) is also described by a n×n q ×n r ×n t complex tensor V a , such that

V a = M m=1 α m V V m ⊗ Q m ⊗ R m ⊗ T m . (41)
Adapting the iterative strategy presented in section 3.4 (the so-called Z-matrix bus method) to the parametric context can be summarized in rewriting [START_REF] Acharya | An analytical approach for DG allocation in primary distribution network[END_REF] with the explicit parametric dependence, i.e.

V [γ+1] a (q, r, t) = Y -1 S * (q, r, t) V * a [γ] (q, r, t) + I 0 . (42)
For algorithmic purposes, and following the ideas presented in section 3.4, this operation is split into two steps: i)First, an intermediate quantity I is computed such that

I(q, r, t) = S * (q, r, t) V * a [γ] (q, r, t) , (43)
then, ii)the second step consists in solving the global (but linear) system, that is in computing

V [γ+1] a (q, r, t) = Y -1 (I(q, r, t) + I 0 ) . ( 44
)
Note that in [START_REF] Frank | Optimal power flow: a bibliographic survey I[END_REF], the operations are not as trivial as in their algebraic version. For instance, for S and

V a in C n , computing I = S * V a is a simple division for each component: [I] l = [S] * l /[V ] * l for l = 1, .
. . , n. For V a (q, r, t) and S(q, r, t) represented in the separable forms [START_REF] Favuzza | Optimal electrical distribution systems reinforcement planning using gas micro turbines by dynamic ant colony search algorithm[END_REF] and [START_REF] Fang | A stochastic power flow method based on polynomial normal transformation and quasi monte carlo simulation[END_REF], the operation (43) requires solving a PGD problem. That is (for each iteration γ) to solve a problem of the type: find I(q, r, t) such that I(q, r, t) V * a [γ] (q, r, t) = S * (q, r, t). The standard PGD procedure consists in computing sequentially the terms of the PGD expansion of I(q, r, t) (loop on M ) and for each term iterate in the alternated directions scheme (this is going to be denoted as a loop on k). The PGD solver uses a greedy algorithm to compute these terms in the expansion (40) (or its tensorial form (41)), see [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | A Short Review on Model Order Reduction Based on Proper Generalized Decomposition[END_REF][START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF][START_REF] Chinesta | Pgd-based computational vademecum for efficient design, optimization and control[END_REF]. The second step is straightforward, since the matrix Y does not depend on the parameters q, r and t, it is possible to conclude that the voltage inherits the same parametric modes of the current, while the vector coefficients are multiplied by Y -1 , for all m = 1, . . . , M .

Thus, in this context, the PGD algorithm involves three nested loops:

-The external one corresponds to the nonlinear solver and iterates in γ -The second is the greedy part of the PGD algorithm to solve (43) (loop on the number of terms of the PGD expansion M ) -The inner loop iterates (for k = 1, 2, . . .) in the alternated direction scheme for each of the parametric dimensions.

The initial solution is typically provided after the slack node intensity, I 0 , namely

V [0] a = Y -1 I 0 . (45)
The global idea of the PDG procedure is very well described and illustrated in [START_REF] García-Blanco | A reduced order modeling approach for optimal allocation of distributed generation in power distribution systems[END_REF] and [START_REF] Borzacchiello | Introduction to the proper generalized decomposition for the solution of the parametrized power equations[END_REF].

Parametric version of the error assessment

As in the case of the algebraic problem, this section is oriented to focus on the methodology to assess the error for the parametric version of the problem.

Although, diverse parameters can be considered, for the sake of simplicity and without losing generality, we present the equations for the parameter r. Note that the behavior of the q and t parametric dimensions is analogous to the r dimension.

Taking the tensorial representation of the solution V a , the error and the residual are also complex matrices in C n×nr ,

E = V -V a , (46)
R(V) = S * -V * (YV -I 0 ). ( 47
)
Now, QoI is taken as the integration with respect to the parametric dimensions of some nonparametric QoI l(•), namely

L(V a (r)) = r l(V a (r)) dr. ( 48
)
Integration respect to the parameter r is determined by the mass matrix M r (associated with the 1D mesh discretizing the parameter r) multiplied by vector 1 nr = [1, 1, . . . , 1] T ∈ C nr . This is possible because the equivalence between the functional and tensorial representations in [START_REF] Favuzza | Optimal electrical distribution systems reinforcement planning using gas micro turbines by dynamic ant colony search algorithm[END_REF] and [START_REF] Florentin | Adaptive reduced basis strategy based on goal oriented error assessment for stochastic problems[END_REF]:

L(V a ) = 1 T nr M r l(V a ), (49)
where l(•) is now the generalization to the tensor introduced in equation [START_REF] Chinesta | Pgd-based computational vademecum for efficient design, optimization and control[END_REF]. The result produces a vector of n r components, namely

l(V a ) = diag((V * a T Y L V a ) Re ), (50)
where the operator diag(•) maps the elements of the diagonal of the input matrix of size n r × n r into a column vector of size n r . Similarly as in [START_REF] De Leon | Iterative solvers in the Newton power flow problem: preconditioners, inexact solutions, and partial jacobian updates[END_REF],

l(V) = l(V a ) + diag((FE) Re ) + diag((GE * ) Re ) ( 51
)
where

F = V * a T Y L and G = (Y L V a )
T are matrices in C nr×n . Using now the Cartesian representation, this equation is rewritten as:

l(V) = l(V a ) + diag( λT p Ê) (52)
where

λp = F Re + G Re -F Im + G Im , Ê = E Re E Im ∈ R 2n×nr .
Using the tensor contraction notation, equation (52) becomes:

L(V) = L(V a ) + λT : Ê ( 53
)
where λ = (1 2n 1 T nr M r ) λp ∈ R 2n×nr . The error equation is derived following the same ideas as in section 3.5:

V a Y E + E * (YV a -I 0 ) = R(V a ). ( 54
)
Taking every column of the matrix V a , it is possible to build two tensors 54) is rewritten as:

A(•, •, ) = Diag(V * a (•, ))Y and B(•, •, ) = Diag((YV a (•, ) -I 0 (•, )))Y for = 1, . . . , n r in C n×n×nr . Thus, (
A • E + B • E * = R(V a ) ( 55
)
where the operation

• denotes a contraction of one index and a Hadamard product in another index. For instance, in the particular case of A ∈ C n×n×nr and E ∈ C n×nr , it reads

A • E i = n j=1
A ij E j , with no sum on . [START_REF] He | An improved power flow algorithm for distri-bution networks based on Z-bus algorithm and forward/backward sweep method[END_REF] Note that the definition is general for any field and for the dimensions of the tensors, the only restriction being that the two last indices of tensor A have the same range as the the two indices of tensor E.

Using the Cartesian representation, the equation becomes linear:

C • Ê = R(V a ), (57)
where

R = E Re R Im ∈ R 2n×nr .
The dual problem is readily introduced as:

C T • ρ = λ, (58)
where

C T (•, •, ) = C(•,
•, ) T , ∀ (transposing only the two first dimensions of the tensor). Hence the error in the quantity of interest using equation ( 57) is:

E QoI = L(V) -L(V a ) = λT : Ê = λ : (C † • R(V a )) = ρT : R(V a ) (59)
where

C † (•, •, ) = C -1 (•,
•, ), ∀ (sectionally inverting the two first dimensions of the tensor).

The error assessment technique using the solution ρT of (58) and the error representation (59) is, in practice, computationally unaffordable. This is due to the multidimensional character of both ρT and R(V a ), which are tensors of order n × n q × n r × n t . Moreover, once ρT and R(V a ) are obtained, all the tensorial dimensions must be contracted (this requires four nested loops) to compute the scalar quantity E QoI .

In the following, we introduce a numerical strategy that condensates all the parametric dimensions in order to devise an amenable error assessment methodology. In this regard, the QoI introduced in ( 48) is integrating the effect of the parametric dimensions and the resulting problem depends only on the physical dimension (represented here by the vector of voltages of size n). Accordingly, the error representation is expected to have the form

E QoI = (ρ A ) T RA (V a ), (60)
where ρA and RA (V a ) are vectors in R 2n that have to be obtained condensing the parametric dimensions (here, integrating with respect to parameter r).

The condensation of R(V a ) ∈ R 2n×nr and C ∈ R 2n×2n×nr into RA (V a ) ∈ R 2n and C A ∈ R 2n×2n (superscript A is used to denote that the quantities are condensed into an accumulated value) is readily obtained by just integrating the parametric dimension, namely

RA (V a ) = r R(•, r) dr = R(V a )M r 1 nr , (61)
and

C A = r C(•, •, r) dr = C M r 1 nr .
It is assumed that there exists some vector ÊA ∈ R 2n , representing an average value of Ê(•, r), such that

r C(•, •, r) • Ê(•, r) dr = C A ÊA . (62)
Consequently, the equation for the mean error ÊA is precisely the following linear system of dimension 2n

C A ÊA = RA (V a ). ( 63
)
Note that the existence of vector ÊA is guaranteed by the integral Mean Value Theorem applied to the lefthand-side of (57), under the hypothesis of having a continuous dependence on Ê(•, r) on r. In this case, there exists some value of r such that ÊA = Ê(•, r). Note that continuity of Ê(•, r) is ensured by the continuity of the parametric description of the solution V a (r). If the modes are not continuous, the existence of ÊA is also guaranteed provided that C A is a regular matrix. In this case, ÊA does not necessarily coincide with any value of Ê(•, r).

In the parametric case, the error in the QoI reads

E QoI =L(V) -L(V a ) = r diag( λp (•, r) T Ê(•, r)) dr = diag( λT p Ê)M r 1 nr , (64)
where the last term in the right uses the multidimensional tensor structure to express the integrals along the r range by a scalar product. An accumulated value of λp , λA ∈ R 2n , is readily introduced λA = r λp (•, r) dr = λp M r 1 nr .

In order to obtain a suitable error representation, it must be assumed that the following hypothesis is true.

Assumption 1 The quantity of interest E QoI is expressed using the accumulated value of λp and the vector ÊA , that is to say,

E QoI = ( λA ) T ÊA .
This can be interpreted as a new application of the mean value theorem in [START_REF] Kabir | Inclusion of slack bus in newton raphson load flow study[END_REF], with the additional assumption that the average value of Ê is again ÊA . Actually, in this case there is not a unique average value: there exists an affine space of dimension 2n -1 where lie all the possible vectors ÊA fulfilling the equation defined in 1. Thus, the assumption claiming that ÊA from equation (62) fulfils also [START_REF] Kabir | Inclusion of slack bus in newton raphson load flow study[END_REF] (at least approximately) is very likely to hold. This assumption is further supported by noting that the dependence on r of C and λp is directly given by the dependence on r of V a (the matrices F, G, A and B and the tensors A and B depend on the solution V a linearly). Thus, the dominant r mode in V a is going to be the dominant r mode also in C and λp and hence ÊA from equation ( 62) is expected to fulfil also [START_REF] Kabir | Inclusion of slack bus in newton raphson load flow study[END_REF]. An error indicator is introduced in section 5.1 in order to numerically check the validity of Assumption 1.

Hence, the dual problem in the condensed form reads

C A T ρA = λA , (65)
and the corresponding error representation is

E QoI = (ρ A ) T RA (V a ) . (66)
Thus, also in the parametric form of the problem, the error in the quantity of interest can be affordably assessed by solving the condensed dual problem [START_REF] Kim | Dispersed Generation Planning Using Improved Hereford Ranch Algorithm[END_REF] and computing the error estimate using [START_REF] Kim | Dispersed generator placement using fuzzyga in distribution systems[END_REF].

The aim of introducing the goal-oriented error estimations in the algorithm presented in section 4.3 is to control the accuracy of the approximation solution through the incorporation of stopping criteria into the procedure. An extended version of the development of the above equations is presented in [START_REF] García-Blanco | Monitoring a PGD solver for parametric power flow problems with goal-oriented error assessment (article in review)[END_REF].

Numerical examples

In this section, the proposed methodology has been applied to the IEEE 57 bus network. The model is onephase grid consisting of 7 generators, 57 buses interconnected by 63 lines and 17 transformers, thus the number of degrees of freedom is n = 57. This network contains PV nodes, although we converted them into PQ nodes without loss of generality by fixing the power vector S. Otherwise, another new loop must be added to the implemented method in order to control the power in the PV nodes, more details about this idea are given in [START_REF] Borzacchiello | Unified formulation of a family of iterative solvers for power systems analysis[END_REF]. The diagram of the network is shown in Figure 1 while the data system can be found in [START_REF]Power systems test case archive[END_REF].

The main objective in the below examples is to solve an optimization problem: find the optimal location and power of a generator that minimizes the system losses, quantity of interest in this work. Since the test system is a transmission one, the concept of generation changes from distributed generator to just generator. As a first step, we compute the solution while we evaluate it for calculating the losses. Secondly, in the post-processing of the results, we look for the value that optimizes the problem statement. In the first step the error assessment is taking into account in the implementation of the solver using some stopping criteria defined in section 5.1. The original data is given in p.u. and we have selected a 100 MVA base for the system structure. It is known that all the buses can be considered as candidates for installing generators and the capacity of the units were chosen between zero and 2 MW.

Tolerances and stopping criteria.

The goal-oriented error estimated defined in the previous section are used to define the stopping criteria in the algorithm described in section 4.2.

Let us introduce the error indicators ξ , where accounts for the type of error measured ( = R for a purely residual estimate; = S for a measure of the stationarity in the loop or = QoI for the error in the quantity of interest), and denotes the loop where it is used ( = γ; = M ; or = k). Thus, the differ-ent stopping criteria are expressed as: continue with the loop while ξ > tol , tol being the different tolerances prescribed for the different criteria.

The definitions of the different error indicators are listed below:

1. Loop in γ ξ R γ = RA (V [γ+1] a ) 2 ŜA 2 , ξ S γ = V [γ+1] a -V [γ] a 2 V [γ+1] a 2
and

ξ QoI γ = |(ρ A ) T RA (V [γ+1] a )| |L(V [γ+1] a )| ( 67
)
where ŜA = S M r 1 nr and • 2 stands for either the L 2 -norm or the Frobenius norm (depending on whether the argument is a vector or a matrix).

Loop in

M ξ R M = RA I (I a ) 2 ŜA 2 , ξ S M = |α M I | |α 1 I | and ξ QoI M = | λT (V [γ],M a -V [γ],M -1 a )| |L(V [γ] a )| , ( 68
)
where RA

I (I a ) = (S -V a I a )M r 1 nr . 3. Loop in k ξ S1 k = (V M ) k+1 -(V M ) k 2 (V M ) k+1 2 and ξ S2 k = (R M ) k+1 -(R M ) k 2 (R M ) k+1 2 . (69)
Moreover, in order to check the stabilization of ρ, the following indicator is introduced:

d ρ = ρ[γ+1] -ρ[γ] 2 ρ[γ+1] 2 . ( 70
)
If the value of d ρ is small enough, the assumption on the stability of ρA is going to be confirmed. Besides, for checking that the Assumption 1 holds, another indicator is introduced:

e Ê = |E QoI -( λA ) T ÊA | |E QoI | . ( 71
)
Note that ÊA is computed using equation ( 63) straightforwardly.

Similarly, the verification of the obtained solution and the corresponding losses is performed with the following error measures (with respect to a reference solution V):

e V = V -V a 2 V 2 (72)
e L = l(V) -l(V a ) 2 l(V) 2 . ( 73
)

Algebraic approach

Here we show a case study where an algebraic version of the power flow problem, associated with the parameters q, r and t, is solved. Particularly, the location of the generator is fixed in the node q = 25, its power is r = 2 p.u. and the hour of the year t = 3454. Consequently, n q , n r , n t = 1 and the number of degrees of freedom is n = 57. Figure 2 illustrates how the error in the quantity of interest ξ QoI γ barely changes when the solution of the dual problem ρ is calculated until the tolerance for the indicator d ρ (in this case 10 -4 ) is reached. Note that we introduce the notation ξ QoI γ for indicating that the vector ρ is computed until the tolerance is reached. The standard notation ξ QoI γ implies that the dual problem is solved at every single iteration. Moreover, Figure 3 shows that the stability of ρ is evident. Since the same fact was noticed in other simulations, from now on in the examples below, once we reached the tolerance for the indicator d ρ , the vector ρ is reused in the following iterations. As a result, at every γ iteration significant computational time is saved. Hence at some point, the cost of calculating the error in the quantity of interest has the same computational cost as the residual calculation because we do not need to update C, λ or ρ.

16

Raquel García-Blanco et al.

Algebraic approach

Here we show a case study where an algebraic version of the power flow problem, associated with the parameters q, r and t, is solved. Particularly, the location of the generator is fixed in the node q = 25, its power is r = 2 p.u. and the hour of the year t = 3454. Consequently, n q , n r , n t = 1 and the number of degrees of freedom is n = 57. Figure 2 illustrates how the error in the quantity of interest ⇠ QoI barely changes when the solution of the dual problem ⇢ is calculated until the tolerance for the indicator d ⇢ (in this case 10 4 ) is reached. Note that we introduce the notation ⇠ QoI ˆ for indicating that the vector ⇢ is computed until the tolerance is reached. The standard notation ⇠ QoI implies that the dual problem is solved at every single iteration. Moreover, Figure 3 shows that the stability of ⇢ is evident. Since the same fact was noticed in other simulations, from now on in the examples below, once we reached the tolerance for the indicator d ⇢ , the vector ⇢ is reused in the following iterations. As a result, at every iteration significant computational time is saved. Hence at some point, the cost of calculating the error in the quantity of interest has the same computational cost as the residual calculation because we do not need to update C, ˆ or ⇢. In order to show the e ciency of the procedure for linearizing the residual and the losses equation explained in section 3.5, the e↵ectivity index is plotted in Figure 4 comparing the relative error with respect to the reference solution e L . Note that for computing the relative errors, that is to say, e L or e V , we consider as real solution the one calculated using Newton-Raphson algorithm while tolerances tol R , tol S and tol QoI are 10 8 .

Parametric approach

Our goal in this section is to solve the optimization problem defined in section 1 of finding the optimal values of the parameters q, r and t when generators are set in the test system. Thus, a parametric version of the power flow problem is involved. In order to show the efficiency of the procedure for linearizing the residual and the losses equation explained in section 3.5, the effectivity index plotted in Figure 4 comparing the relative error with respect to the reference solution e L . Note that for computing the relative errors, that is to say, e L or e V , we consider as 16

Raquel García-Blanco et al.

Algebraic approach

Here we show a case study where an algebraic version of the power flow problem, associated with the parameters q, r and t, is solved. Particularly, the location of the generator is fixed in the node q = 25, its power is r = 2 p.u. and the hour of the year t = 3454. Consequently, n q , n r , n t = 1 and the number of degrees of freedom is n = 57. Figure 2 illustrates how the error in the quantity of interest ⇠ QoI barely changes when the solution of the dual problem ⇢ is calculated until the tolerance for the indicator d ⇢ (in this case 10 4 ) is reached. Note that we introduce the notation ⇠ QoI ˆ for indicating that the vector ⇢ is computed until the tolerance is reached. The standard notation ⇠ QoI implies that the dual problem is solved at every single iteration. Moreover, Figure 3 shows that the stability of ⇢ is evident. Since the same fact was noticed in other simulations, from now on in the examples below, once we reached the tolerance for the indicator d ⇢ , the vector ⇢ is reused in the following iterations. As a result, at every iteration significant computational time is saved. Hence at some point, the cost of calculating the error in the quantity of interest has the same computational cost as the residual calculation because we do not need to update C, ˆ or ⇢. In order to show the e ciency of the procedure for linearizing the residual and the losses equation explained in section 3.5, the e↵ectivity index is plotted in Figure 4 comparing the relative error with respect to the reference solution e L . Note that for computing the relative errors, that is to say, e L or e V , we consider as real solution the one calculated using Newton-Raphson algorithm while tolerances tol R , tol S and tol QoI are 10 8 .

Parametric approach

Our goal in this section is to solve the optimization problem defined in section 1 of finding the optimal values of the parameters q, r and t when generators are set in the test system. Thus, a parametric version of the power flow problem is involved. real solution the one calculated using Newton-Raphson algorithm while tolerances tol R γ , tol S γ and tol QoI γ are 10 .

16

Raquel García-Blanco et al.

Algebraic approach

Here we show a case study where an algebraic version of the power flow problem, associated with the parameters q, r and t, is solved. Particularly, the location of the generator is fixed in the node q = 25, its power is r = 2 p.u. and the hour of the year t = 3454. Consequently, n q , n r , n t = 1 and the number of degrees of freedom is n = 57. Figure 2 illustrates how the error in the quantity of interest ⇠ QoI barely changes when the solution of the dual problem ⇢ is calculated until the tolerance for the indicator d ⇢ (in this case 10 4 ) is reached. Note that we introduce the notation ⇠ QoI ˆ for indicating that the vector ⇢ is computed until the tolerance is reached. The standard notation ⇠ QoI implies that the dual problem is solved at every single iteration. Moreover, Figure 3 shows that the stability of ⇢ is evident. Since the same fact was noticed in other simulations, from now on in the examples below, once we reached the tolerance for the indicator d ⇢ , the vector ⇢ is reused in the following iterations. As a result, at every iteration significant computational time is saved. Hence at some point, the cost of calculating the error in the quantity of interest has the same computational cost as the residual calculation because we do not need to update C, ˆ or ⇢. In order to show the e ciency of the procedure for linearizing the residual and the losses equation explained in section 3.5, the e↵ectivity index is plotted in Figure 4 comparing the relative error with respect to the reference solution e L . Note that for computing the relative errors, that is to say, e L or e V , we consider as real solution the one calculated using Newton-Raphson algorithm while tolerances tol R , tol S and tol QoI are 10 8 .

Parametric approach

Our goal in this section is to solve the optimization problem defined in section 1 of finding the optimal values of the parameters q, r and t when generators are set in the test system. Thus, a parametric version of the power flow problem is involved.

Parametric approach

Our goal in this section is to solve the optimization problem defined in section 1 of finding the optimal values of the parameters q, r and t when generators are set in the test system. Thus, a parametric version of the power flow problem is involved.

Optimal nominal power and location of a generator with fixed loads.

As a first example, the parameter t is fixed while q and r vary. Hence, we seek the value of the parameters q and r that minimize the system losses. Since the loads are fixed, the problem consists in finding the voltage solution as a separated representation:

V a = M m α m V V m ⊗ Q m ⊗ R m , (74)
where r belongs to a set of possible values of power that the generator can provide. That is, the partition of the interval [0, r max ] where the increment is r max /(n r -1) with r max = 2 p.u and n r is the number of samples, particularly in this example n r = 100. We set the generator in different nodes along the network corresponding with q = 2, . . . , 26, being n q = 25 while n t = 1. Note that the number of degrees of freedom is n × n q × n r = 57 × 25 × 100.

The separated representation of the input data S is:

S = α 1 S S 1 ⊗ 1 nr + α 2 S Q 2 ⊗ r 2 (75)
where S 1 is the vector of demand loads with dimension n, 1 nr is a vector of ones with dimension n r , Q 2 is a zero vector except for the location of the generator where 1 is placed, r 2 is a vector in R nr where [r 2 ] h = (h/n r )•r max , for all h = 0, . . . , n r and α 1 S , α 2 S are positive scalars. The novice of the presented technique is that it makes possible to control the quality of the solution in term of the losses, the quantity of interest in this case, during the iterative process. Thus, the goal is to set different tolerances and compare the obtained solutions in order to validate the goal-oriented error estimates. The first list of tolerances is tol γ = 10 -5 , tol M = 10 -6 and tol k = 10 -7 for = S, R, QoI. Figure 5 and6 show the stopping criteria for γ and M respectively. The numbers in Figure 5 represent the amount of modes that the solution contains at every iteration γ. In this case, the final solution consists of 10 modes after 20 iterations.

The same quantities are shown in Figures 7 and8 but the fixed tolerances are tol γ = 10 -7 , tol M = 10 -9 and tol k = 10 -10 for = S, R, QoI in this case. As we can see, the amount of terms of the solution changes, we need more than the double of terms (a total of 21) in this simulation. In both cases, after a few iterations the number of modes that the solutions contain are stabilized.

The relative errors for the losses e L and the reference solution e V are shown in Figures 9 and 10 respectively. It is observed that the solution that contains more terms is more accurate. This is because at every itle Suppressed Due to Excessive Length 17

.3.1 Optimal nominal power and location of a enerator with fixed loads.

s a first example, the parameter t is fixed while q and vary. Hence, we seek the value of the parameters q nd r that minimize the system losses. Since the loads re fixed, the problem consists in finding the voltage olution as a separated representation:

a = M X m ↵ m V V m ⌦ Q m ⌦ R m , (74)
here r belongs to a set of possible values of power that he generator can provide. That is, the partition of the nterval [0, r max ] where the increment is r max /(n r 1) ith r max = 2 p.u and n r is the number of samples, articularly in this example n r = 100. We set the genertor in di↵erent nodes along the network corresponding ith q = 2, . . . , 26, being n q = 25 while n t = 1. Note hat the number of degrees of freedom is n ⇥ n q ⇥ n r = 7 ⇥ 25 ⇥ 100. The separated representation of the input data S is:

S = ↵ 1 S S 1 ⌦ 1 n r + ↵ 2 S Q 2 ⌦ r 2 (75)
here S 1 is the vector of demand loads with dimension , 1 n r is a vector of ones with dimension n r , Q 2 is a zero ector except for the location of the generator where 1 is laced, r 2 is a vector in R n r where [r 2 ] h = (h/n r )•r max , or all h = 0, . . . , n r and ↵ 1 S , ↵ 2 S are positive scalars. The novice of the presented technique is that it akes possible to control the quality of the solution in erm of the losses, the quantity of interest in this case, uring the iterative process. Thus, the goal is to set diferent tolerances and compare the obtained solutions in rder to validate the goal-oriented error estimates. The rst list of tolerances is tol ⇤ = 10 5 , tol ⇤ M = 10 6 and ol ⇤ k = 10 7 for ⇤ = S, R, QoI. Figure 5 and6 show the topping criteria for and M respectively. The numers in Figure 5 represent the amount of modes the olution contains at every iteration . In this case, the nal solution consists of 10 modes after 20 iterations.

The same quantities are shown in Figures 7 and8 ut the fixed tolerances are tol ⇤ = 10 7 , tol ⇤ M = 10 9 nd tol ⇤ k = 10 10 for ⇤ = S, R, QoI in this case. As we an see, the amount of terms of the solution changes, e need more than the double of terms (a total of 21) n this simulation. In both cases, after a few iterations he number of modes that the solutions contain are stailized.

The relative errors for the losses e L and the refernce solution e V are shown in Figures 9 and 10 respecively. It is observed that the solution that contains ore terms is more accurate. This is because at every iteration M , we add a new term hence more information is considered. Clearly, this added information is enough for changing significantly the quality of the solution.

Based on the numerical results, the introduction of error estimators in the procedure allows to control the whole procedure and specifically the construction of the solution V. But, a priori we neglected some terms in the equations assuming that the achieved result will be accurate enough and also taking into account the above-mentioned Assumption 1. itle Suppressed Due to Excessive Length 17

3.1 Optimal nominal power and location of a nerator with fixed loads.

s a first example, the parameter t is fixed while q and vary. Hence, we seek the value of the parameters q d r that minimize the system losses. Since the loads e fixed, the problem consists in finding the voltage lution as a separated representation:

a = M X m ↵ m V V m ⌦ Q m ⌦ R m , (74)
here r belongs to a set of possible values of power that e generator can provide. That is, the partition of the terval [0, r max ] where the increment is r max /(n r 1) ith r max = 2 p.u and n r is the number of samples, articularly in this example n r = 100. We set the generor in di↵erent nodes along the network corresponding ith q = 2, . . . , 26, being n q = 25 while n t = 1. Note at the number of degrees of freedom is n ⇥ n q ⇥ n r = ⇥ 25 ⇥ 100. The separated representation of the input data S is:

= ↵ 1 S S 1 ⌦ 1 n r + ↵ 2 S Q 2 ⌦ r 2 (75)
here S 1 is the vector of demand loads with dimension , 1 n r is a vector of ones with dimension n r , Q 2 is a zero ctor except for the location of the generator where 1 is laced, r 2 is a vector in R n r where [r 2 ] h = (h/n r )•r max , r all h = 0, . . . , n r and ↵ 1 S , ↵ 2 S are positive scalars. The novice of the presented technique is that it akes possible to control the quality of the solution in rm of the losses, the quantity of interest in this case, uring the iterative process. Thus, the goal is to set difrent tolerances and compare the obtained solutions in der to validate the goal-oriented error estimates. The rst list of tolerances is tol ⇤ = 10 5 , tol ⇤ M = 10 6 and l ⇤ k = 10 7 for ⇤ = S, R, QoI. Figure 5 and6 show the opping criteria for and M respectively. The numers in Figure 5 represent the amount of modes that the lution contains at every iteration . In this case, the nal solution consists of 10 modes after 20 iterations.

The same quantities are shown in Figures 7 and8 ut the fixed tolerances are tol ⇤ = 10 7 , tol ⇤ M = 10 9 d tol ⇤ k = 10 10 for ⇤ = S, R, QoI in this case. As we n see, the amount of terms of the solution changes, e need more than the double of terms (a total of [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF] this simulation. In both cases, after a few iterations e number of modes that the solutions contain are stailized.

The relative errors for the losses e L and the referce solution e V are shown in Figures 9 and 10 respecvely. It is observed that the solution that contains ore terms is more accurate. This is because at every iteration M , we add a new term hence more information is considered. Clearly, this added information is enough for changing significantly the quality of the solution.

Based on the numerical results, the introduction of error estimators in the procedure allows to control the whole procedure and specifically the construction of the solution V. But, a priori we neglected some terms in the equations assuming that the achieved result will be accurate enough and also taking into account the above-mentioned Assumption 1. iteration M , we add a new term hence more information is considered. Clearly, this added information enough for changing significantly the quality of the solution.

Based on the numerical results, the introduction of error estimators in the procedure allows to control the whole procedure and specifically the construction of the solution V. But, a priori we neglected some terms in the equations assuming that the achieved result will be accurate enough and also taking into account the above-mentioned Assumption 1.

18

Raquel García-Blanco et al. Figures 11 and 12 illustrate the validity of these two assumptions based on this numerical example. Specifically, Figure 11 proves that the Assumption 1 holds in the numerical example using the error indicator e Ê while Figure 12 shows the e↵ectivity index comparing the relative error in losses e L and the error indicator of the quantity of interest ⇠ QoI .

Once we get the approximation V a , the losses are calculated using the operator L. The two dimensional representation of the losses is presented in Figure 13 losses corresponds to a generator situated in the node q = 12. The minimum loss is 0.201 p.u. given when the power of the generator is r = 0.78 p.u. The optimization step could be carried out using any algorithm, however since the objective function is now explicitly available, it just required a simple evaluation.

5.3.2

Optimal location of a generator with time varying loads.

In this second example, a parametric power flow problem is solved for seeking the values of q and r that Figures 11 and 12 illustrate the validity of these two assumptions based on this numerical example. Specifically, Figure 11 proves that the Assumption 1 holds in the numerical example using the error indicator e Ê while Figure 12 shows the e↵ectivity index comparing the relative error in losses e L and the error indicator of the quantity of interest ⇠ QoI .

Once we get the approximation V a , the losses are calculated using the operator L. The two dimensional representation of the losses is presented in Figure 13 where we can observe that the minimum value of the losses corresponds to a generator situated in the node q = 12. The minimum loss is 0.201 p.u. given when the power of the generator is r = 0.78 p.u. The optimization step could be carried out using any algorithm, however since the objective function is now explicitly available, it just required a simple evaluation.

5.3.2

Optimal location of a generator with time varying loads.

In this second example, a parametric power flow problem is solved for seeking the values of q and r that minimize losses when the parameter time t is also con- Figures 11 and 12 illustrate the validity of these two assumptions based on this numerical example. Specifically, Figure 11 proves that the Assumption 1 holds in the numerical example using the error indicator e Ê while Figure 12 shows the effectivity index comparing the relative error in losses e L and the error indicator of the quantity of interest ξ QoI γ . Once we get the approximation V a , the losses are calculated using the operator L. The two dimensional representation of the losses is Figures 11 and 12 illustrate the validity of these two assumptions based on this numerical example. Specifically, Figure 11 proves that the Assumption 1 holds in the numerical example using the error indicator e Ê while Figure 12 shows the e↵ectivity index comparing the relative error in losses e L and the error indicator of the quantity of interest ⇠ QoI .

Once we get the approximation V a , the losses are calculated using the operator L. The two dimensional representation of the losses is presented in Figure 13 where we can observe that the minimum value of the losses corresponds to a generator situated in the node q = 12. The minimum loss is 0.201 p.u. given when the power of the generator is r = 0.78 p.u. The optimization step could be carried out using any algorithm, however since the objective function is now explicitly available, it just required a simple evaluation.

5.3.2

Optimal location of a generator with time varying loads.

In this second example, a parametric power flow problem is solved for seeking the values of q and r that minimize losses when the parameter time t is also con- Figures 11 and 12 illustrate the validity of these two ssumptions based on this numerical example. Specifcally, Figure 11 proves that the Assumption 1 holds n the numerical example using the error indicator e Ê hile Figure 12 shows the e↵ectivity index comparing he relative error in losses e L and the error indicator of he quantity of interest ⇠ QoI . Once we get the approximation V a , the losses are alculated using the operator L. The two dimensional epresentation of the losses is presented in Figure 13 here we can observe that the minimum value of the losses corresponds to a generator situated in the node q = 12. The minimum loss is 0.201 p.u. given when the power of the generator is r = 0.78 p.u. The optimization step could be carried out using any algorithm, however since the objective function is now explicitly available, it just required a simple evaluation.

5.3.2

Optimal location of a generator with time varying loads.

In this second example, a parametric power flow problem is solved for seeking the values of q and r that minimize losses when the parameter time t is also con- presented in Figure 13 where we can observe that the minimum value of the losses corresponds to a generator situated in the node q = 12. The minimum loss is 0.201 p.u. given when the power of the generator is r = 0.78 p.u. The optimization step could be carried out using any algorithm, however since the objective function is now explicitly available, it just required a simple evaluation. Fig. 13 Reconstructed System Losses sidered. The solution in this case reads as:

V a = M X m ↵ m V V m ⌦ Q m ⌦ R m ⌦ T m , (76)
while the representation of the data load term S requires more terms:

S = 5 X h=1 ↵ h S Sh ⌦ Qh ⌦ Ȓh ⌦ T h + ↵ 1 S 1 ⌦ Q 1 ⌦ R 1 ⌦ T 1 (77)
where S 1 and Sh are n vectors representing the nodal positions of the network, Q 1 and Qh are zero vectors of n q components except for the location of the generator varying from 25 to 50 in the second term, thus n q = 25, and R 1 and Ȓh corresponds to the variation of the power, for all h = 1, . . . , 5. In this case, the maximum power is r max = 1.5 p.u. and the number of samples is n r = 100. In order to represent the time in the test system, the load demand and the generation profiles during a year are represented by the load curves T h (t), 8h = 1, . . . , 5 for the nodes and T 1 fot the generator. These curves were generated using the software HOMER described in [START_REF] Lambert | Micropower System Modeling with HOMER[END_REF] for the test system described in a former work, see [START_REF] García-Blanco | Monitoring a PGD solver for parametric power flow problems with goal-oriented error assessment (article in review)[END_REF]. The curves vary from 1 to 8760 with a time step of 1h, thus n t = 8760. Following the same procedure as before, we fix the same two shorts list of tolerances. In the first case, the tolerances are tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 and in the second case, these are tol ⇤ = 10 7 , tol ⇤ M = 10 8 and tol ⇤ k = 10 10 for ⇤ = S, R, QoI. The number of degrees of freedom in both cases is n ⇥ n q ⇥ n r ⇥ n t = 57 ⇥ 25 ⇥ 100 ⇥ 8760. That is the reason of why the number of iterations for reaching the tolerances in both examples are higher comparing to the numerical example above, 33 versus 50 iterations when the fixed tolerances are smaller. Accordingly, the number of modes of the final solution also varies, 14 versus 31 modes. This fact is shown in Figures 14,15, 16 and 17 where the convergence diagram of the stopping criteria in the greedy algorithm are also plotted. It is worth mentioning that in all M iterations, at some point the criterion ⇠ R M stabilizes after some iterations. This might be because at every iteration M , we add a new term hence more information is considered. However, it is possible that the added information is not enough for changing significantly the quality of the solution, thus the residual in the first step of the algorithm does not decrease.

The diagrams of the relative errors comparing both approximations with the reference solution and their

V a = M X m ↵ m V V m ⌦ Q m ⌦ R m ⌦ T m , (76)
while the representation of the data load term S requires more terms:

S = 5 X h=1 ↵ h S Sh ⌦ Qh ⌦ Ȓh ⌦ T h + ↵ 1 S 1 ⌦ Q 1 ⌦ R 1 ⌦ T 1 (77)
where S 1 and Sh are n vectors representing the nodal positions of the network, Q 1 and Qh are zero vectors of n q components except for the location of the generator varying from 25 to 50 in the second term, thus n q = 25, and R 1 and Ȓh corresponds to the variation of the power, for all h = 1, . . . , 5. In this case, the maximum power is r max = 1.5 p.u. and the number of samples is n r = 100. In order to represent the time in the test system, the load demand and the generation profiles during a year are represented by the load curves T h (t), 8h = 1, . . . , 5 for the nodes and T 1 fot the generator. These curves were generated using the software HOMER described in [START_REF] Lambert | Micropower System Modeling with HOMER[END_REF] for the test system described in a former work, see [START_REF] García-Blanco | Monitoring a PGD solver for parametric power flow problems with goal-oriented error assessment (article in review)[END_REF]. The curves vary from 1 to 8760 with a time step of 1h, thus n t = 8760. Following the same procedure as before, we fix the same two shorts list of tolerances. In the first case, the tolerances are tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 and in the second case, these are tol ⇤ = 10 7 , tol ⇤ M = 10 8 and tol ⇤ k = 10 10 for ⇤ = S, R, QoI. The number of degrees of freedom in both cases is n ⇥ n q ⇥ n r ⇥ n t = 57 ⇥ 25 ⇥ 100 ⇥ 8760. That is the reason of why the number of iterations for reaching the tolerances in both examples are higher comparing to the numerical example above, 33 versus 50 iterations when the fixed tolerances are smaller. Accordingly, the number of modes of the final solution also varies, 14 versus 31 modes. This fact is shown in Figures 14,15, 16 and 17 where the convergence diagram of the stopping criteria in the greedy algorithm are also plotted. It is worth mentioning that in all M iterations, at some point the criterion ⇠ R M stabilizes after some iterations. This might be because at every iteration M , we add a new term hence more information is considered. However, it is possible that the added information is not enough for changing significantly the quality of the solution, thus the residual in the first step of the algorithm does not decrease.

The diagrams of the relative errors comparing both approximations with the reference solution and their

Optimal location of a generator with time varying loads.

In this second example, a parametric power flow problem is solved for seeking the values of q and r that minimize losses when the parameter time t is also considered. The solution in this case reads as:

V a = M m α m V V m ⊗ Q m ⊗ R m ⊗ T m , (76)
S = 5 h=1 α h S Sh ⊗ Qh ⊗ Ȓh ⊗ T h + α 1 S 1 ⊗ Q 1 ⊗ R 1 ⊗ T 1 (77)
where S 1 and Sh are n vectors representing the nodal positions of the network, Q 1 and Qh are zero vectors of n q components except for the location of the generator varying from 25 to 50 in the second term, thus n q = 25, and R 1 and Ȓh corresponds to the variation of the power, for all h = 1, . . . , 5. In this case, the maximum power is r max = 1.5 p.u. and the number of samples is n r = 100. In order to represent the time in the test system, the load demand and the generation profiles during a year are represented by the load curves T h (t), ∀h = 1, . . . , 5 for the nodes and T 1 fot the generator. These curves were generated using the software HOMER described in [START_REF] Lambert | Micropower System Modeling with HOMER[END_REF] for the test system described in a former work, see [START_REF] García-Blanco | Monitoring a PGD solver for parametric power flow problems with goal-oriented error assessment (article in review)[END_REF]. The curves vary from 1 to 8760 with a time step of 1h, thus n t = 8760.

Following the same procedure as before, we fix the same two shorts list of tolerances. In the first case, the tolerances are tol γ = 10 -5 , tol M = 10 -6 and tol k = 10 -7 and in the second case, these are tol γ = 10 -7 , tol M = 10 -8 and tol k = 10 -10 for = S, R, QoI. The number of degrees of freedom in both cases is n × n q × n r × n t = 57 × 25 × 100 × 8760. That is the reason of why the number of γ iterations for reaching the tolerances in both examples are higher comparing to the numerical example above, 33 versus 50 iterations when the fixed tolerances are smaller. Accordingly, the number of modes of the final solution also varies, 14 versus 31 modes. This fact is shown in Figures 14,15 It is worth mentioning that in all M iterations, at some point the criterion ξ R M stabilizes after some iterations. This might be because at every iteration M , we add a new term hence more information is considered. However, it is possible that the added information is not enough for changing significantly the quality of the solution, thus the residual in the first step of the algorithm does not decrease. that the solution that contains more modes is signifi- The diagrams of the relative errors comparing both approximations with the reference solution and their losses are shown in Figures 18 and19. It is observed that the solution that contains more modes is significantly more accurate than the one obtained when the tolerances are higher.

When the parameter t is involved in the parametric power flow problem, the losses are annual, that means, at each node the operador L is the sum of the losses at every hour of the year. Hence, the reconstruction of the losses shown in Figure 20 is two dimensional. The location of the generator that provides a minimal losses, 0.712 p.u., is q = 23 while the power is r = 1.2 p.u.

Conclusions

In this paper we presented a review of methods involved in solving the power flow problem in both algebraic and parametric versions. Additionally, a brief summary of our contributions in this field is also introduced. that the solution that contains more modes is significantly more accurate than the one obtained when the tolerances are higher.

When the parameter t is involved in the parametric power flow problem, the losses are annual, that means, at each node the operador L is the sum of the losses at every hour of the year. Hence, the reconstruction of the losses shown in Figure 20 is two dimensional. The location of the generator that provides a minimal losses, 0.712 p.u., is q = 23 while the power is r = 1.2 p.u.

Co

In th in so para our T braic but genc volta The classification of solvers for the power flow algebraic equation is based on their historical appearance, but also on their main characteristics: rate of convergence, robustness (in terms of converging to the high voltage solution or the computational time) and memory storage during the procedure. The application of the Alternating Search Directions methods applied to this version of the problem provides a family of methods able to recover some classical ones as NR or GS.

Conclusions

In this paper we presented a review of methods involved in solving the power flow problem in both algebraic and parametric versions. Additionally, a brief summary of our contributions in this field is also introduced. The classification of solvers for the power flow algebraic equation is based on their historical appearance, but also on their main characteristics: rate of convergence, robustness (in terms of converging to the high voltage solution or the computational time) and mem- In both versions of the problem, a goal-oriented error estimation strategy is introduced. This technique provides consistent stopping criteria for all the iterative schemes of the PGD algorithm. The error is based on the residual and the losses equation which have to be linearized neglecting the quadratic terms. Furthermore, a Cartesian representation is needed since the conjugate function is not holomorphic.

The standard strategy in the context of the error assessment requires the definition of an adjoint (or dual) problem. In practice, as it can be seen in the numerical examples, it is observed to be stationary along the iterative process saving computational resources. Another Besides, it achieves a good performance in terms of accuracy and computational time. Parametric solvers are analyzed, mainly the PLF and OPF, pointing out specifically their competency for solving optimization problems when generators are set in the networks. In order to avoid the curse of dimensionality when Parametric Power Flow problem is solved, Reduced Order Models are introduced. Particularly, PGD technique is the strategy adopted for solving the parametric version of the problem since it presents a double advantage: the separated representation of the In both versions of the problem, a goal-oriented error estimation strategy is introduced. This technique provides consistent stopping criteria for all the iterative schemes of the PGD algorithm. The error is based on the residual and the losses equation which have to be linearized neglecting the quadratic terms. Furthermore, a Cartesian representation is needed since the conjugate function is not holomorphic.

The standard strategy in the context of the error assessment requires the definition of an adjoint (or dual) problem. In practice, as it can be seen in the numerical examples, it is observed to be stationary along the iterative process saving computational resources. Another singular characteristic of the proposed error strategy is that the parametric dimensions need to be condensed. This is done integrating respect to the corresponding parameters the error and the losses equation. It is essential to assume that accumulated error values coincide after integrating both equations. The numerical examples show that this Assumption is accurately fulfilled.

This new computational strategy is used in optimization problems when generators are set in a grid, based on losses minimization. The novelty of this procedure is that the admissible error in the losses is fixed a priori, therefore the solution is built adding the necessary terms in its PGD representation.

Fig. 1

 1 Fig. 1 Single-line diagram of standard IEEE 57 bus test system

Fig. 2

 2 Fig. 2 Convergence diagram of the error in the quantity of interest with the iteration index .

Fig. 3

 3 Fig. 3 Convergence diagram of stagnation criteria for the solution of the dual problem ⇢.

Fig. 4

 4 Fig. 4 E↵ectivity index in the losses.

Fig. 2

 2 Fig. 2 Convergence diagram of the error in the quantity of interest with the iteration index γ.

Fig. 2

 2 Fig. 2 Convergence diagram of the error in the quantity of interest with the iteration index .

Fig. 3

 3 Fig. 3 Convergence diagram of stagnation criteria for the solution of the dual problem ⇢.

Fig. 4

 4 Fig. 4 E↵ectivity index in the losses.

Fig. 3

 3 Fig. 3 Convergence diagram of stagnation criteria for the solution of dual problem ρ.

Fig. 2

 2 Fig. 2 Convergence diagram of the error in the quantity of interest with the iteration index .

Fig. 3

 3 Fig. 3 Convergence diagram of stagnation criteria for the solution of the dual problem ⇢.

Fig. 4

 4 Fig. 4 E↵ectivity index in the losses.

Fig. 4

 4 Fig. 4 Effectivity index in the losses.

Fig. 5

 5 Fig. 5 Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 6

 6 Fig. 6 Convergence diagram of the enrichment algorithm in the last M iteration (last -iteration) for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 .

Fig. 5

 5 Fig. 5 Convergence diagram of the stopping criteria for outer loop with the index γ for tolerances tol γ = 10 -5 , tol M = 10 -6 and tol k = 10 -7 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 5 Fig. 6

 56 Fig. 5 Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 6

 6 Fig. 6 Convergence diagram of the enrichment algorithm in the last M iteration (last γ-iteration) for tolerances tol γ = 10 -5 , tol M = 10 -6 and tol = 10 -7 .

Fig. 7 Fig. 8

 78 Fig. 7 Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 7 , tol ⇤ M = 10 9 and tol ⇤ k = 10 10 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 9 Fig. 10

 910 Fig. 9 Relative errors comparing the real and the approximated solution for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 .

Fig. 7 Fig. 7 Fig. 8

 778 Fig. 7 Convergence diagram of the stopping criteria for the outer loop with the iteration index γ for tolerances tol γ = 10 -7 , tol M = 10 -9 and tol k = 10 -10 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 9 Fig. 10

 910 Fig. 9 Relative errors comparing the real and the approximated solution for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 .

Fig. 8

 8 Fig. 8 Convergence diagram of the enrichment algorithm in the last M iteration (last γ-iteration) for tolerances tol γ = 10 -7 , tol M = 10 -9 and tol k = 10 -10 .

Fig. 7 Fig. 8

 78 Fig. 7 Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 7 , tol ⇤ M = 10 9 and tol ⇤ k = 10 10 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 9 Fig. 10

 910 Fig. 9 Relative errors comparing the real and the approximated solution for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 .

Fig. 9 7 8

 978 Fig. 9 Relative errors comparing the real and the approximated solution for tolerances tol γ = 10 -5 , tol M = 10 -6 and tol k = 10 -7 .

Fig. 9 Fig. 10

 910 Fig. 9 Relative errors comparing the real and the approximated solution for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 .

Fig. 10

 10 Fig. 10 Relative errors comparing the real and the approximated solution for tolerances tol γ = 10 -7 , tol M = 10 -9 and tol k = 10 -10 .

Fig. 11 ⇠

 11 Fig. 11 Convergence diagram of the error indicator e Ê with the iteration index .

Fig. 11 Fig. 11 ⇠Fig. 13

 111113 Fig. 11 Convergence diagram of the error indicator e Ê with the iteration index γ.

Fig. 12

 12 Fig. 12 Effectivity index in the losses.

Fig. 13

 13 Fig. 13 Reconstructed System Losses

  , 16 and 17 where the convergence diagram of the stopping criteria in the greedy algorithm are also plotted.

Fig. 14 1 MFig. 15

 14115 Fig. 14 Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 16

 16 Fig. 16 Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 7 , tol ⇤ M = 10 9 and tol ⇤ k = 10 10 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 17 6 ConclusionsFig. 14

 17614 Fig. 17 Convergence diagram of the enrichment algorithm in the last M iteration (last -iteration) for tolerances tol ⇤ = 10 7 , tol ⇤ M = 10 9 and tol ⇤ k = 10 10 .

Fig. 14 Fig. 15

 1415 Fig.[START_REF] Brameller | Some improved methods for digital network analysis[END_REF] Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 . The numbers along the curves refer to the number of modes that the solution contains.

  Fig. ou 10 th co

Fig. 15

 15 Fig. 15 Convergence diagram of the enrichment algorithm in the M iteration (last γ-iteration) for tolerances tol γ = 10 -5 , tol M = 10 -6 and tol k = 10 -7 . 20 Raquel García-Blanco et al.

Fig. 14

 14 Fig. 14 Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 16

 16 Fig. 16 Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 7 , tol ⇤ M = 10 9 and tol ⇤ k = 10 10 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 16

 16 Fig.[START_REF] Brown | Z-matrix algorithms in load-flow programs[END_REF] Convergence diagram of the stopping criteria for the outer loop with the iteration index γ for tolerances tol γ = 10 -7 , tol M = 10 -9 and tol k = 10 -10 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 16 Fig. 17

 1617 Fig.[START_REF] Brown | Z-matrix algorithms in load-flow programs[END_REF] Convergence diagram of the stopping criteria for the outer loop with the iteration index for tolerances tol ⇤ = 10 7 , tol ⇤ M = 10 9 and tol ⇤ k = 10 10 . The numbers along the curves refer to the number of modes that the solution contains.

Fig. 17 Fig. 18 Fig. 20

 171820 Fig. 17 Convergence diagram of the enrichment algorithm in the last M iteration (last γ-iteration) for tolerances tol γ = 10 -7 , tol M = 10 -9 and tol k = 10 -10 . Title Suppressed Due to Excessive Length 21

Fig. 18

 18 Fig. 18 Relative errors comparing the real the approximated solution for tolerances tol γ = 10 -5 , tol M = 10 -6 and k = 10 -7 .

Fig. 18 Fig. 19 Fig. 19 1 GenFig. 20

 181919120 Fig. 18 Relative errors comparing the real and the approximated solution for tolerances tol ⇤ = 10 5 , tol ⇤ M = 10 6 and tol ⇤ k = 10 7 .

Table 1

 1 Search directions for the classical methods

Acknowledgements The first two authors acknowledge the financial support of Ministerio de Economía y Competitividad, grant DPI2014-51844-C2-2-R.

On behalf of all authors, the corresponding author states that there is no conflict of interest.