Observability Analysis and Restoration for Distribution Grids: A Robust Greedy Numerical Algorithm

This paper presents a robust and efficient algorithm for analyzing and restoring the observability of distribution grids. A linearized load flow model is used to characterize three- and multi-phase distribution grids under balanced and unbalanced conditions. Assuming that metering devices are installed to measure nodal magnitudes at some buses, the rank of the matrix associated with the non-metered magnitudes is determined by using a greedy strategy combined with a Gram-Schmidt orthonormalization procedure. Unlike other numerical-based techniques, the proposed algorithm avoids common errors caused by denominators close to zero. As an additional salient feature, a mathematically sound stopping criterion is provided to guarantee the robustness of the proposed algorithm. Furthermore, for unobservable grids, the algorithm identifies the linearly dependent rows in the aforementioned matrix, which correspond to the buses that should be metered to restore grid observability. The effective performance of the proposed algorithm is illustrated using several medium- and low-voltage benchmarks with up to 906 buses.


Observability Analysis and Restoration for Distribution Grids: A Robust Greedy Numerical Algorithm
Natalia Alguacil , Senior Member, IEEE, Henar Herrero , and Cristina Solares Abstract-This paper presents a robust and efficient algorithm for analyzing and restoring the observability of distribution grids.A linearized load flow model is used to characterize threeand multi-phase distribution grids under balanced and unbalanced conditions.Assuming that metering devices are installed to measure nodal magnitudes at some buses, the rank of the matrix associated with the non-metered magnitudes is determined by using a greedy strategy combined with a Gram-Schmidt orthonormalization procedure.Unlike other numerical-based techniques, the proposed algorithm avoids common errors caused by denominators close to zero.As an additional salient feature, a mathematically sound stopping criterion is provided to guarantee the robustness of the proposed algorithm.Furthermore, for unobservable grids, the algorithm identifies the linearly dependent rows in the aforementioned matrix, which correspond to the buses that should be metered to restore grid observability.The effective performance of the proposed algorithm is illustrated using several medium-and low-voltage benchmarks with up to 906 buses.Index Terms-Distribution grids, multi-phase networks, observability analysis, observability restoration, phasor measurement unit (PMU).

I. INTRODUCTION
O BSERVABILITY analysis and restoration have become essential for the reliable and cost-efficient operation of the increasingly active distribution grids [1].In the past, due to their passive behavior and the scarcity of metering devices, these networks were typically operated with limited observability [2], [3].However, with the recent deployment of smart meters [4] and phasor measurement units (PMUs) [5], full observability of distribution grids can be achieved or restored Natalia Alguacil is with the Departamento de Ingeniería Eléctrica, Electrónica, Automática y Comunicaciones, E.T.S.I.Industrial, Universidad de Castilla-La Mancha, 13071 Ciudad Real, Spain (e-mail: Natalia.Alguacil@ uclm.es).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TSG.2023.3256488.
Digital Object Identifier 10.1109/TSG.2023.3256488if such devices are conveniently located [6].Furthermore, the observability of distribution grids can be improved if the complementary strengths of smart meters and PMUs are exploited [7].
Similarly to the observability analysis carried out for transmission systems, accurate and updated information on their topology, network parameters, and metering locations is also required in distribution systems [3], [8].However, since distribution grids are reconfigured more frequently than transmission networks [9], computing limitations are tighter for these grids.As a consequence, dynamic observability analysis should be efficiently implemented to capture topological changes during grid operation.In addition to the reconfiguration issue, distribution grids may operate under unbalanced conditions caused by both untransposed distribution lines and unbalanced loads, which requires a multi-phase modeling [10].In such cases, the phases in which single-or multi-phase equipment is connected should be clearly identified [11].Note that phase identification is part of topology identification [12], [13], which in turn is required for observability analysis.
In [23], grid observability was studied using a topological algorithm in which the radial nature of the network was exploited.Likewise, a topological algorithm was proposed in [24], in which the rank of the system matrix is determined by analyzing the topology of the system graph.A numerical criterion was also provided to analyze the rank deficiency of the system matrix.
The numerical algorithm described in [25] was based on the inverse function theory and the Jacobian matrix of the system.Moreover, an orthogonal linear rotational transformation was used to address the large r/x ratio in distribution systems.In [26], a multi-phase numerical observability analysis was developed using an orthogonal factorization of the three-phase measurement Jacobian matrix.To that end, a weak coupling P-v and Q-θ , as in transmission systems, was considered.The work presented in [27] proposed a numerical algorithm for observability analysis based on the transformation of the Gain matrix using Cholesky decomposition.A decoupled PQ grid model typically used in transmission systems was revisited to identify unobservable branches and critical data for unbalanced three-phase distribution systems.
In addition, a triangular factorization of the Gain matrix was performed considering branch current magnitudes and angles as the state variables.In [28], a methodology to generate pseudo-measurements along the distribution network was proposed to guarantee full observability of the grid.Recently, in [29], correlations between measurements were taken into account to improve system observability and an observability metric useful for state estimation was provided.Note that in many studies the effect of the distribution network is modeled by linear load flow equations such as those presented in [30].It is also worth emphasizing that the computational burden of existing techniques [23], [24], [25], [26], [27], [28], [29] grows with the size of the required vectors and matrices, which may be up to three times larger for unbalanced operation as compared to a balanced setting.Thus, one of the yet unresolved challenges for both topological and numerical algorithms is related to their scalability when addressing large distribution grids.Moreover, observability restoration is typically neglected in the related literature, [27] being a relevant exception.Unfortunately, the approach described in [27] is prone to issues associated with its reliance on the method presented in [14], which may involve dividing by numbers close to zero and/or pivoting positions.
As an alternative, a novel robust greedy algorithm is proposed here for both observability analysis and observability restoration of three-and multi-phase distribution grids under balanced and unbalanced conditions.Based on a customarily used linear network model [30], the proposed algorithm starts by conveniently partitioning the system matrix, as done in [18], in which the problem of observability in linear systems of equations was addressed.Bearing in mind that measurements of active/reactive power injections, voltage magnitudes, and voltage angles are only available for a subset of buses, the algorithm determines whether these measurements are sufficient or not to guarantee grid observability.That is, if the corresponding nodal magnitudes in the non-metered buses can be inferred from the available measurements.Furthermore, in case of unobservability, the algorithm identifies the buses at which meters should be located to restore observability while considering practical multi-phase metering devices.Several case studies, including multi-phase grids [31], are examined to comprehensively illustrate the performance in small instances and the computational efficiency of the algorithm in larger instances.
This paper significantly departs from the existing literature in the following methodological aspects.First, the proposed algorithm calculates the rank of the system matrix, an optimal basis of the row space, and the coordinates of the remaining rows in that basis.Interestingly, the methodology requires neither inverting matrices nor selecting pivots, while avoiding the errors caused when dividing by values close to zero [25], [26], [27], [28], [29].Moreover, unlike the algorithms considered in [25], [26], [27], [28], [29], the proposed approach is numerically stable since it does not require setting any tolerance to distinguish between zero and nonzero values.Additionally, in contrast to [23], [24], [25], [26], [28], and [29] observability restoration is considered.Furthermore, the proposed approach differs from the algorithm presented in [22] for transmission grids in two crucial aspects: 1) the algorithm relies on a mathematically sound stopping criterion that does not require tuning tolerances, and 2) the restoration of observability is carried out by identifying the best candidate buses to be metered using a practical criterion suitable for distribution grids.Finally, the above distinctive features allow analyzing and restoring observability for instances that are far larger than those addressed in [23], [24], [25], [26], [27], [28], [29].
The main contributions of this work are as follows: • Based on an existing numerical algorithm for observability analysis in transmission systems and a widely used linear load flow model, a novel observability approach is developed to accommodate the singularities of both balanced and unbalanced distribution grids.As a major novelty within the context of distribution grids, the proposed methodology avoids the errors caused by dividing by values close to zero.• Unlike its transmission-related counterpart and most approaches for distribution systems, the restoration of grid observability is handled by identifying the nonmetered buses at which smart meters or PMUs should be located.Furthermore, consistently with practical metering devices that meter all the phases of a bus, observability is restored by including the best candidate bus in the subset of metered buses with all its phases as a single block.• A mathematically sound stopping criterion is provided to guarantee the robustness of the algorithm, thus avoiding tuning tolerances.• The efficiency of the algorithm is empirically proven using two IEEE benchmarks and a 906-bus feeder under different scenarios, which are generated by modifying both the number and the location of metered buses.Thus, for the first time in the literature, distribution system operators are provided with a practical tool to analyze and restore observability.The rest of this paper is organized as follows.In Section II, the observability matrices for distribution systems are derived.In Section III, the proposed algorithm for observability analysis and restoration is presented.In Section IV, numerical results from several case studies are provided and discussed.In Section V, the main conclusions are summarized.Finally, a theorem is included in the Appendix to justify the suitability of the stopping criterion of the proposed algorithm.

II. PROBLEM STATEMENT
This section presents the matrices required for the proposed algorithm for distribution grid observability analysis and restoration.Based on [24] and [30], the following linear approximation of the load flow equations is considered: where v is the vector of voltage magnitudes, θ is the vector of voltage angles, p is the vector of active power nodal Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
injections, and q is the vector of reactive power nodal injections, whereas, for a given grid topology, R and X are the real-valued matrices computed from the line resistances and reactances, respectively, as done in [24].

A. Single-Phase Formulation With Smart Meter Data
Under balanced conditions, three-phase distribution grids can be characterized using an equivalent single-phase model.Thus, the vectors and matrices in ( 1) and ( 2) are n-and n × ndimensional, respectively.When smart meters are deployed, the nodal vectors v, p, and q are available.As done in [24], it is assumed that smart meters are located at a subset of buses M ⊂ N , where N = {1, 2, . . ., n} is the set of all buses in the distribution grid.Thus, the above nodal vectors can be partitioned as and q = [q M , q M], where M = N \M denotes the subset of non-metered buses.
Using the notation introduced above, the system of equations in (1) can be partitioned as follows: where subindices M and M refer to the rows/columns of the matrices and vectors associated with the metered and nonmetered buses, respectively.System (3) can be rearranged to place the non-metered magnitudes in the left-hand side and the metered magnitudes in the right-hand side.That is, where 0 denotes a matrix of zeroes and I is the identity matrix.
To check if the system is observable, i.e., if v, p, and q can be determined for the non-metered buses, the transpose of the matrix in the left-hand side of (4) is required [22], which is formulated as follows: Matrix H is analyzed using the iterative procedure described in Section III.

B. Single-Phase Formulation With PMU Data
If PMUs are installed in the grid, expressions ( 1) and ( 2) should be jointly considered since v, p, q, and θ are metered by the available PMUs.To that end, both systems of equations are partitioned as described above.Thus, the partitioning of ( 1) yields (3) whereas the system of equations ( 2) is partitioned analogously giving rise to: Next, equations ( 3) and ( 6) are conveniently reordered to obtain the following system of equations: Similar to (5), matrix H can be obtained by transposing the matrix in the left-hand side of ( 7):

C. Multi-Phase Formulation
The H matrices in ( 5) and ( 8) corresponding to the equivalent single-phase formulations can be readily extended to the general multi-phase setting.To that end, the dimensions of the vectors and matrices for unbalanced three-phase systems should be increased as follows: v is the 3n-dimensional vector of voltage magnitudes, θ is the 3n-dimensional vector of voltage angles, p is the 3n-dimensional vector of active power nodal injections, and q is the 3n-dimensional vector of reactive power nodal injections.Analogously, R and X are real-valued 3n × 3n-dimensional matrices [30].Note that, for multi-phase systems with components with either 1, 2, or 3 phases, the dimensions of the above vectors and matrices should be adjusted accordingly.These modifications pave the way for the development of a novel algorithm to analyze and restore the observability of multi-phase distribution grids.

III. PROPOSED ALGORITHM
As sketched in the pseudocode below, the proposed algorithm, which is referred to as Algorithm 1, involves two phases, namely, Analysis Phase, and Restoration Phase.Note that the main steps of the algorithm are labeled at the left hand-side of the pseudo-codes.Both phases are described as in Algorithms 2 and 3.
In the Analysis Phase the maximum number of row vectors in matrix H that are linearly independent are selected, so that a basis of the subspace of I R l generated by the selected row vectors of H is obtained.The selection process requires orthonormalizing the basis to calculate distances as explained next.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

A4
Orthonormalize the vectors {w 1 , h s 2 } by Gram-Schmidt and obtain vector w 2 ; An arbitrary row vector in H, h k , is first selected and normalized as w 1 = h k / h k .The next vector selected is the farthest from the subspace generated by vector w 1 , denoted by W 1 .To identify which is the farthest vector, the distance of a generic row vector h j to the subspace W 1 is calculated as

R3
Add the bus corresponding to z max and all its phases to the subset of metered buses M;

R4
Update matrix H and the corresponding values of m and l accordingly; return H; That is, the distance of vector h j to subspace W 1 requires the projection of this vector on the subspace, which, in turn, is defined through the orthonormal basis of W 1 .Once the farthest vector is found, the subspace, W 2 , generated by this newly identified vector and vector w 1 is considered.Similarly, the farthest vector from this new subspace is selected by first orthonormalizing the basis formed by vectors w 1 and w 2 , and then projecting the other row vectors on W 2 so that the corresponding distances to this subspace W 2 are calculated.This process is repeated by selecting new vectors.Note that the maximum distance decreases iteratively.The process stops when either the matrix is full rank or the maximum distance is greater than that of the previous iteration.For a theoretical justification of this stopping criterion, the interested reader is referred to Theorem 1 in the Appendix.
By considering the farthest row vector at each iteration, the row dependencies are efficiently identified.If the distance of the farthest row vector to a subspace is 0, this vector is linearly dependent on the vectors already included in the subspace and so are the remaining row vectors outside the subspace.Note that a distance equal to 0 can only be attained once all linearly independent vectors have been identified.As backed by Theorem 1, this situation occurs when the maximum distance at one iteration is greater than that of the previous iteration.Since this stopping criterion does not require setting an arbitrary tolerance to distinguish between zero and nonzero values, we claim that the proposed procedure is numerically stable as opposed to existing ones.
The output of the Analysis Phase is: • The rank of matrix H, i.e., N.
• Vector s, whose entries are integers indexing the rows in H that are linearly independent.• The basis S N , which is a set formed by the orthonormal vectors obtained throughout the iterative process.
, whose columns are the coordinates of the rows of H expressed in the basis S N .If N = m, i.e., if the number of components in s is equal to m, the system is observable.Otherwise, if N < m, matrix A is the starting point to restore grid observability.
As for the proposed Restoration Phase, unlike [22], where row dependencies were not exploited, observability is restored by analyzing the columns of matrix C, where is a submatrix of A. The entries of vector t are the indexes of the rows in H that are linearly dependent, so that these entries identify the columns in A that are selected to be included in matrix C. Note that each row in C is associated with a vector in the set S N while each column corresponds to a linearly dependent row in H expressed in the basis S N .In particular, the original criterion proposed in this work to decide which is the best candidate bus to be metered consists in finding the column vector of C with the largest number of coordinates in S N less than .Once the algorithm identifies a bus as the best candidate, all its phases are included in the subset of metered buses M. In case that none of the column vectors of C satisfies the -criterion, any bus associated with the columns of C can be arbitrarily selected to be included in M. It is worth emphasizing that the selection of the bus to be added to M takes into consideration that practical metering devices in distribution networks are multi-phase.It is remarkable that observability is successfully restored regardless of the selected value for , which only impacts on the selection of buses, and, hence, on the number of iterations in the Restoration Phase.
The main differences of the analysis procedure with respect to the algorithm presented in [22] lie in Step A10 and in the Restoration Phase.As per Step A10, the Analysis Phase terminates when the proposed stopping criterion e I ≥ e I−1 is satisfied.Note that this condition was not used in [22].It should be also emphasized that this work improves upon [22] by providing the aforementioned theoretical justification for this novel criterion, namely, Theorem 1 in the Appendix.
It should be highlighted that the proposed algorithm is related to a QR factorization, but the rows of matrix H are rearranged in a different order by using an optimal orthonormalization procedure.The proposed technique is based on the Gram-Schmidt orthonormalization of a selected basis using a greedy strategy over the space generated by the rows of matrix H.The algorithm is greedy as it makes the best choice at each step [32].Note that rows are not sequentially chosen as done under the classical Gram-Schmidt orthonormalization.Rather, at each iteration, the row with the largest distance to the corresponding subspace is selected.Interestingly, despite its greedy nature, no approximation is involved.
Finally, it should be noted that the computational complexity of the proposed algorithm is O(N 3 ).

IV. CASE STUDIES
To illustrate the effective performance of the proposed algorithm for analyzing and restoring grid observability, several benchmarks are examined under balanced and unbalanced conditions, assuming that multi-phase metering devices are available for the metered buses.First, simulations for the IEEE 4-bus feeder [33] are run for didactical purposes.Then, observability is checked and restored for the IEEE 13-bus feeder [34], which is a medium-voltage multi-phase distribution grid.Finally, the computational efficiency of the algorithm is illustrated on the European 906-bus low-voltage feeder [10].For all instances, is set to 10 −7 according to the order of magnitude of the involved numerical values.Simulations were run using MATLAB [35] on a workstation SIE Ladon with two Intel Xeon Gold 6248R processors at 3 GHz and 768 GB of RAM.

A. IEEE 4-Bus Feeder
The first case study considers a four-wire line configuration and a substation transformer rated 12.47 kV/4.16kV and 6 MVA.For expository purposes, different operating conditions and metering scenarios are used to analyze the observability of this grid.
1) Balanced Operation: Balanced conditions are examined first.Thus, the single-phase (1φ) formulations presented in Sections II-A and II-B are used.Fig. 1 shows the one-line diagram for a representative metering scenario wherein smart meters are located at buses 1 and 3.Under this scenario, which is referred to as scenario 1φ-a, the magnitudes p 1 , q 1 , v 1 , p 3 , q 3 , and v 3 are metered and the corresponding matrix H is presented in (9).Note that labels p 2 , q 2 , and v 2 in (9) denote the correspondence between the non-metered magnitudes and the rows of this matrix.
Considering the rows of matrix H in (9), the Analysis Phase is applied as follows.In Step A1, the first row of H is arbitrarily selected.Once this row is normalized, vector w 1 =[0.2143,0.6907, 0.6907] is obtained.In Step A2, the projections of rows 2 and 3 on the space W 1 = span{w 1 } are calculated.Then, the distances of each projection are computed as E (1)
It should be emphasized that, for this benchmark, nonobservability is identified for any single-phase scenario in which bus 1 is non-metered, no matter whether the remaining buses are metered or not.Fig. 2 shows a different scenario, namely, scenario 1φ-b, in which the magnitudes p 2 , q 2 , v 2 , p 3 , q 3 , and v 3 are non-metered.For this scenario, matrix H is given in (10).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
At the third iteration of the Analysis Phase (I = 3), the largest distance is equal to e 3 = 1.415 × 10 −16 .Moreover, three linearly independent vectors are obtained and, hence, matrix H has no full rank.The output of the Analysis Phase at this iteration is: • The rank of H, N, which is equal to 3, i.e., N = l.
• Vector s = (1, 5, 6), whose entries are the indices of the rows of matrix H that are linearly independent.• The basis S 3 given in (11), whose rows are the orthonormal vectors that are obtained along the algorithm.
Once Step R1 is applied, vector t = (2, 3, 4) is obtained thus yielding matrix C in (12).This matrix includes the coordinates of the linearly dependent rows in H expressed in the basis S 3 .
Thus, under scenario 1φ-b, the grid is deemed as nonobservable since matrix H in ( 10) is rank-deficient.
To restore observability, matrix C in ( 12) is explored.Each vector in S 3 , which is a row in the matrix given in (11), is associated with a linearly independent row in H.Then, as denoted by the labels in C, the rows in this matrix correspond to the linearly independent magnitudes in (10), i.e., p 2 , v 2 , and v 3 , whereas the columns in C are related to the linearly dependent magnitudes in H, namely, p 3 , q 2 , and q 3 .In this case, the -criterion in Step R2 is not met since all the values in C are larger than 10 −7 .That is, it is indifferent to select either bus 2 or bus 3. Let us consider that a smart meter can be placed at bus 2. Adding bus 2 to the subset of metered buses results in a new matrix H with all rows linearly independent, thereby identifying the system as observable.Analogously, if a smart meter is placed at bus 3, the system would also be deemed as observable.
This benchmark is also useful to illustrate the single-phase analysis when PMUs are available.To that end, scenario 1φ-c is considered, wherein buses 1 and 3 are non-metered.For this scenario, the proposed algorithm is initially fed with the matrix H given in (13).As a result, matrix H is identified as rankdeficient and six linearly independent vectors are obtained.0.0000 − 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 − 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 − 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 − 1.0000 To restore observability, the application of Step R1 of the algorithm yields the matrix C given in (14).
According to the -criterion, the algorithm selects the second column of matrix C in ( 14) since it features the maximum number of entries with an absolute value less than 10 −7 , namely one.Note that this column is associated with the reactive power injected at bus 1.Thus, bus 1 is added to the subset of metered buses and matrix H is updated accordingly.Since the new matrix H is of full rank, observability can be restored if a PMU is placed at bus 1.
2) Unbalanced Operation: The use of the multi-phase formulation presented in Section II-C is illustrated with a scenario considering unbalanced operation.Under this scenario, it is assumed that no smart meter is placed at bus 2, i.e., the magnitudes p 2 , q 2 , and v 2 are not metered in all phases.
For this case, matrix H is given in (15), shown at the bottom of the next page.The application of the proposed approach to this multi-phase instance gives rise to the following outcome: • The rank of H is 9.
• Vector s = (1, 9, • The basis S 9 in ( 16), shown at the bottom of the next page.Thus, the IEEE 4-bus feeder for this particular multiphase scenario is deemed as observable since the unknown magnitudes p 2 , q 2 , and v 2 in all phases can be uniquely determined.

B. IEEE 13-Bus Feeder
The IEEE 13-bus feeder has laterals with one, two, and three phases.The delta-wye substation transformer is rated 115 kV/4.16kV and 5 kVA.
Considering balanced operation and smart meters, the algorithm has been first applied to all possible instances wherein a single bus is non-metered.As done in [36], the single-phase equivalent feeder is obtained by assuming that distribution lines are symmetrically transposed and every bus is supplied by all three phases.From this study, it can be concluded that Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the IEEE 13-bus feeder is unobservable when buses 1, 5, and 8 are the only non-metered buses.
Next, the metering scenario 1φ-a depicted in Fig. 3 is evaluated.For this scenario, it is assumed that magnitudes p i , q i , and v i are known for all buses except buses 2 and 7.The Analysis Phase requires six iterations to find six linearly independent rows in the corresponding matrix H, with the largest distance being e 6 = 1.0281 × 10 −2 .Thus, the grid is identified as observable under this scenario since matrix H is of full rank (N = m = 6), and, hence, p 2 , q 2 , v 2 , p 7 , q 7 , and v 7 can be uniquely determined.
The proposed approach has also been applied to the metering scenario 1φ-b shown in Fig. 4, with p i , q i , and v i , for i = 2, 3, 7, 12, being non-metered.For this scenario, the stopping criterion is satisfied for I = 11, with e 11 = 2.9302 × 10 −5 > e 10 = 1.669 × 10 −12 and N = 10.Thus, the grid is deemed as non-observable since only 10 linearly independent rows are found in H, which is a rank-deficient matrix.
To restore observability, matrix C in ( 17) is computed.Note that the rows/columns in C are associated with the linearly independent/dependent rows in H.
The second column of C in (17), which is associated with the active power injected at bus 7, features the largest number of zeroes.Thus, bus 7 is added to the subset of metered buses and matrix H is updated accordingly.For this new H, the grid is still identified as non-observable.As a consequence, bus 3 is added to the subset of metered buses in the next iteration to eventually restore grid observability.0.0353 0.2487 q 3 0.0853 −0.0049 q 12 −0.0021−0.0302 q 2 −0.0125 0.0000 p 12 0.0000 0.0000 (17) Moreover, a new scenario is generated assuming that PMUs are available at all buses except for buses 2, 3, 7, and 12, for which the corresponding magnitudes p i , q i , v i , and θ i are unknown.After 16 iterations of the analysis procedure, the proposed approach determines that this matrix has full rank with p i , q i , v i , and θ i being uniquely determined for i = 2, 3, 7, 12. Therefore, the grid is observable for this scenario.

TABLE I IEEE 13-BUS FEEDER -RESULTS
Additional representative metering scenarios with either smart meters or PMUs have been specifically designed to further illustrate the performance of the proposed approach.In Table I, the results for the single-and multi-phase models are summarized.It should be noted that the computing times reported in this table correspond to the observability analysis, as observability restoration involves the installation of extra devices, which is to be decided in the future.
As evidenced in Table I, the grid is observable for scenarios in which the percentage of non-metered buses(-phases) ranges between 11% and 33%.Moreover, the percentage of non-metered buses(-phases) that should be metered to restore observability ranges between 20% and 60%.Computing times are similar in all scenarios, ranging between O(10 −3 ) s and O(10 −1 ) s.
Finally, in order to validate the proposed methodology, all metering scenarios considered in [24] have been analyzed under balanced conditions and identical results are obtained, as expected.

C. European 906-Bus Low-Voltage Feeder
In the European 906-bus low-voltage feeder (ELV), all laterals have three phases and the substation transformer, rated 11 kV/416 V and 800 kVA, has a delta-wye configuration.For this benchmark, both single-and multi-phase models are analyzed.For the sake of simplicity, the metering scenarios are randomly generated including only smart meters or PMUs.
Table II lists the size of matrix H and the required computing time per iteration for metering scenarios with full observability.Note that the computational burden grows with the matrix size.III summarizes the results for metering scenarios requiring observability restoration.It should be noted that scenarios 1-4 are associated with the single-phase formulation while scenario 5 corresponds to the multi-phase model.
In scenario 1, with 45 non-metered buses (5% of the total number of buses), the algorithm finds that 26 rows of H are linearly dependent and that 25 buses, i.e., 56% of the non-metered buses, are required to restore observability.In scenario 2, with the same number of non-metered buses as scenario 1 but with PMUs, the algorithm finds that only 2 rows of H are linearly dependent and that one bus (50% of the nonmetered buses) is required to restore observability.In scenario 3, for which the number of non-metered buses is increased to 181 (20% of the total number of buses), the algorithm finds that 166 rows of H are linearly dependent and that 147 buses, i.e., 81% of the non-metered buses, are required to restore observability.In scenario 4, with the same number of nonmetered buses as scenario 3 but with PMUs, the algorithm finds that 26 rows of H are linearly dependent and that 11 buses (6% of the non-metered buses) are required to restore observability.Finally, in the multi-phase scenario 5 with 91 non-metered buses (10% of the total number of buses), the From the results summarized in Table III, it can be concluded that the number of non-metered buses required to restore observability is scenario-dependent.It is also worth mentioning that the percentage of metered buses is not correlated with the percentage of buses in which metering devices are required to restore observability.Additionally, it should be emphasized that, for the scenarios under consideration, the size of the matrices H ranges between 135 × 905 and 819 × 2715, while the corresponding computing times are between 4 s and 9005 s.
Finally, the application of the stopping criterion in Step A10 is illustrated using scenario 1 in Table III.This instance, which is not observable, has 109 linearly independent rows in H.The maximum distances e I at iterations I = 108, 109, and 110 are, respectively, 4 × 10 −4 , 2 × 10 −11 , and 6 × 10 −4 .It should be noted that e 109 is not exactly zero due to rounding errors.As can be noted, the maximum distance decreases down to a value near zero at iteration I = 109, whereas it increases at iteration I = 110.Additionally, according to Theorem 1 in the Appendix, pr S 109 h j = h j , j = 1, . . ., 135, which in turn implies that E (109) j = ||h j − pr S 109 h j || = 0 and h j ∈ W 109 , j = 1, . . ., 135.Thus, the maximum number of linearly independent rows is N = 109.In other words, as the maximum distance increases at iteration I = 110, the 110 th vector is linearly dependent on previous vectors.Thus, at iteration I = 110 the Analysis Phase of the algorithm terminates with N = 109 providing the information required to start the Restoration Phase.
The stopping criterion proposed in [22], i.e., e I ≤ 10 −7 , has been alternatively considered, yielding identical results for this scenario.However, if a different tolerance is the algorithm may require a different number of iterations and even provide an incorrect output.In particular, for e I ≤ 10 −15 , the Analysis Phase terminates for N = 135 with matrix H being of full rank, thereby identifying the grid as observable.This behavior is depicted in Fig. 5, in which the logarithm with base 10 of the largest distance e I is plotted for each iteration I. From this figure, it is inferred that the vectors between iterations 110 and 135 are included in the basis, thus requiring 25 additional iterations as compared to the proposed stopping criterion and leading to a wrong observability analysis.
Therefore, checking if the maximum distance decreases down to a minimum value and then stops decreasing at the next iteration is an effective stopping criterion.This finding is theoretically backed by Theorem 1 in the Appendix.

V. CONCLUSION
This paper has presented a robust greedy algorithm for analyzing and restoring observability in distribution grids.Unlike previous works, the algorithm proposed to determine the rank of the system matrix is numerically stable.The observability analysis phase relies on the iterative determination of the number of linearly independent rows in the matrix representing system operation.If needed, observability restoration is subsequently incorporated as a relevant distinctive feature while using an effective practical criterion for selecting the best candidate buses to be metered.Moreover, the robustness of the algorithm is guaranteed by a mathematically sound stopping criterion that does not require setting any tolerance to distinguish between zero and nonzero values.Thus, the proposed numerical algorithm overcomes the limitations of existing approaches by avoiding divisions by values close to zero.
Furthermore, the proposed algorithm is suitable for threephase and multi-phase distribution grids under both balanced and unbalanced conditions.Results from several case studies based on medium-and low-voltage distribution grids with up to 906 buses demonstrate the effective performance and the practical applicability of the proposed algorithm.
Further research will consider the generalization of the proposed formulation to accommodate other metering devices.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Manuscript received 10
June 2022; revised 18 October 2022 and 30 January 2023; accepted 8 March 2023.Date of publication 14 March 2023; date of current version 23 October 2023.This work was supported in part by Grant PID2021-126566OB-I00 and Grant PID2021-122579OB-I00, funded by the Spanish Ministry of Science and Innovation MCIN/AEI/10.13039/501100011033and by the ERDF, EU "A Way of Making Europe"; and in part by Grant SBPLY/21/180501/000154, funded by the Junta de Comunidades de Castilla-La Mancha, Spain, and by the ERDF.Paper no. TSG-00822-2022.(Corresponding author: Natalia Alguacil.)

Fig. 5 .
Fig.5.ELV -Evolution of log 10 of the maximum distance along the iterative process using an alternative stopping criterion.