An efficient method to include equality constraints in branch current distribution system state estimation

Distribution system state estimation is a fundamental tool for the management and control functions envisaged for future distribution grids. The design of accurate and efficient algorithms is essential to provide estimates compliant with the needed accuracy requirements and to allow the real-time operation of the different applications. To achieve such requirements, peculiarities of the distribution systems have to be duly taken into account. Branch current-based estimators are an efficient solution for performing state estimation in radial or weakly meshed networks. In this paper, a simple technique, which exploits the particular formulation of the branch current estimators, is proposed to deal with zero injection and mesh constraints. Tests performed on an unbalanced IEEE 123-bus network show the capability of the proposed method to further improve efficiency performance of branch current estimators.


Introduction
In the smart grid (SG) scenario, where control and management activities of the electric distribution network are expected to play a relevant and increasing role, distribution system state estimation (DSSE) is conceived as a fundamental monitoring tool. In fact, control systems, such as distribution management systems (DMSs), must rely on a possibly complete and accurate knowledge of the state of the network given by DSSE [1].
In the current evolving scenario, the increasing presence of distributed generation (DG), storage devices and flexible loads to be controlled leads to unforeseen dynamics, which require suitable measurement responsiveness. In this context, DSSE techniques able to work at high reporting rates, while safeguarding the estimation accuracies required by specific applications, are needed. DSSE has to be able to include all the different measurement types, provided by traditional and modern measurement devices, which can be available with different frequencies and accuracies. In particular, the phasor measurement *Correspondence: carlo@diee.unica.it Department of Electrical and Electronic Engineering, University of Cagliari, Piazza d'Armi 1, 09123 Cagliari, Italy units (PMUs) [2,3], which give phasor measurements of both voltages and currents synchronized with respect to a common time reference (the so-called synchrophasors) are becoming increasingly widespread in transmission systems [4][5][6] and are expected to be widely used also in distribution systems, along with new-generation power meters.
For these reasons, innovative and dedicated solutions able to estimate the operating point of the future distribution network are increasingly needed.
DSSE techniques are based on measurement models that link measurements with the state variables to be estimated. They are intended to elaborate the measurements acquired from the field, by means of distributed measurement systems, along with all the a priori information that can be collected about load and generator activity.
In fact, it is impractical and economically infeasible to have a fully monitored distribution network, where each node is equipped with a measurement device connected to the monitoring infrastructure. The DSSE thus reaches observability by relying on the so-called pseudomeasurements [7], which include historical or forecast data on generator production and load consumption. DSSE approaches proposed in the literature are mainly based on weighted least squares (WLS) algorithms [7][8][9][10][11][12] and mainly differ between each other in the chosen state variables and in the way the different measurement types are included. In particular, two main classes of DSSE exist, based on two different choices of the underlying state vector: node voltage state estimators, NV-DSSEs (as in [7,13]) and branch current estimators, BC-DSSEs (as [8][9][10][11][12]). BC-DSSE is suitably designed to better keep into account the peculiar characteristics of distribution systems, as the radial or weakly meshed topology and the high r/x ratios, and it is usually faster with respect to those based on voltage state (see [14]).
An important topic, which can influence both accuracy and speed of DSSE based on WLS methods, concerns the constraints given by a priori knowledge on network operation and topology. For instance, a priori information exploitable by DSSE includes also the identification of the so-called 'zero injections', that is of the nodes that are surely known to have no power consumption or generation. Zero injections are frequent in a distribution grid, in particular, because in a three-phase unbalanced context, some nodes may have no loads or generators connected to some of the phases. Besides, the inactivity of a load or generator, if it represents an absolutely sure information, could be also translated in a zero injection constraint. Additional constraints can be also present. As an example, in a branch current formulation, possible meshes have to be duly considered.
Such constraints can be treated in different ways. The simplest method is to consider them as virtual pseudomeasurements and to use a large weight to enforce them. This choice can lead to numerical conditioning problems [15], and thus other options have been advanced in the literature. In [16], for instance, it is proposed to include them using constraints expressed through Lagrange multipliers, while in [17], they are treated as normal measurements with a low weight and the constraints are re-imposed between the subsequent iterations of the WLS.
In this paper, a simple way to deal with the equality constraints, well-suited to BC-DSSE (and in particular to the efficient formulation presented in [12]) and based on state vector reduction, is proposed. This approach is compared with other traditional and commonly used techniques to underline the advantages by means of simulation results obtained on a IEEE 123-buses three-phase test network.

Branch current state estimation
State estimation techniques are based on mathematical relations between system state variables and measurements collected from the distributed measurement system. The measurements in a distribution grid can be the traditional ones, as voltage and current magnitudes, real and reactive power flows, and power injections at buses, or the current and voltage synchrophasors provided by PMUs. Usually, distribution networks are only partially monitored. As a consequence, prior information on the loads (the so-called pseudo-measurements) are necessary to perform state estimation. Thus, the forecasts of power injections usually constitute the majority of the measurements available for DSSE.
As aforementioned, different state variables can be considered for the estimation algorithm and in particular node voltages or branch currents. Such variables can be represented either in polar or rectangular coordinates. In this paper, the enhanced branch current-based estimator proposed in [12] is adopted, as it was shown to be as accurate as those traditionally based on node voltages and more efficient in the practical application of distribution networks. This algorithm will be referred to in the following as BC-DSSE. The general measurement model adopted for state estimation is: The state vector x, in the BC-DSSE, is given by the branch currents of all the N br network branches, in rectangular coordinates, and the voltage v s at a reference node, for instance the slack bus. In a three-phase framework, , under the hypothesis that synchronized measurements are present. This formulation exploits the absolute phase angles provided by PMU measurements, which use the common time reference of the coordinated universal time (UTC). Such time synchronization can be obtained by means of Global Positioning system (GPS) or other synchronization sources (see for instance [18,19]).
Pseudo-measurements, in the model (1), are handled as measurements that are assigned with a higher standard deviation σ to highlight the lower accuracy due to the fact they are not based on real measurements but rather on historical and forecast data.
In BC-DSSE, the estimation of the state x is obtained by an iterative algorithm. Each iteration consists of: • definition/update of measurements and residuals; • branch current estimation applying a WLS method; • network voltage state computation through a forward sweep calculation.
For each iteration k, in the first step of BC-DSSE, power measurements are translated into equivalent phasor current measurements using the node voltages estimated in the previous iteration. This approach allows including power measurements (and above all pseudomeasurements) easily in the estimator, since equivalent current injections are linearly linked to the branch current variables. Using the updated vector z k of the measurements and the previously estimated statex k−1 , the measurement residuals r k = z k − h(x k−1 ) are computed. In the first iteration, when estimates are not still available, an initialization of the state variables is needed.
The WLS step is then used, at each iteration, to find the state variable variation x k =x k −x k−1 that minimizes the weighted sum of the squares of the residuals, by solving the following normal equations: where H k = H(x k−1 ) is the Jacobian of the measurement functions at iteration k and W is the weighting matrix, equal to the inverse of z . Matrix G = H T WH (subscript k will be dropped in the following, for the sake of simplicity) represents the socalled gain matrix, which has to be inverted or factorized to find the solution of (2). Such matrix and its characteristics play thus a key role in the estimation process.
In the last step of the algorithm, network voltages for each node are computed, by a simple evaluation of voltage drops along the lines, starting from the estimated voltage in the reference bus. With a matrix expression, the node voltage phasors v k , at iteration k, are obtained as follows: where Z paths is the matrix that contains, for each row i, the branch impedances z j that belong to the path that links v i to the reference bus v s . It is worth noting that, in the case of meshed networks, a radial tree of the network can be considered in order to identify the paths linking the slack bus to each node of the grid and to make the forward sweep step possible. The chosen tree can be whatever, since the inclusion of the mesh constraints in the preceding WLS step ensures the final achievement of the same voltage results independently from the particular choice of the path.
The procedure is repeated until a given threshold in estimated state variation is reached.

Formulation of the equality constraints
Zero injections are the most common case of equality constraints that can be found in a distribution system. They can be generally represented by: where c(·) is a N c size vector that represents the constraints to be kept into account in state estimation. It is worth noting that these constraints can be nonlinear, depending on the chosen state variables of the system, but they are linear in the case of rectangular BC-DSSE. In a similar way, possible presence of meshes also leads to equality constraints that, in the branch current formulation, have to be suitably considered. The way to handle such constraints in the DSSE can affect the performance of the estimator, in terms of both accuracy and speed. Different methods have been used in the literature. In the following, virtual measurements and Lagrange multipliers are first described, in order to present the equality constraints issue in a self-contained discussion. Then, a new method, based on state vector reduction, is proposed.

Virtual measurement method
Zero injections can be treated as virtual measurements given by (4), that is as measurements to be included into the measurement model (1) with a very low measurement uncertainty (represented by the virtual standard deviation σ zi ). Such approach leads to a WLS where a high weight is attributed to the power measurement of zero injection nodes. It is important to remark that the resulting weight is not representative of a real uncertainty (which should be equal to zero), but it is adopted only to point out the higher reliability of the information associated to the virtual measurements with respect to the other available measurements. In the DSSE, the measurement vector becomes z = z T m , z T zi T , where z m and z zi are the proper measurement vector and the virtual measurement vector, respectively. Indicating with C the Jacobian of the zero injection constraints, and using the previous notation for the weighting matrix and the residuals, they can also be divided as follows: The zero injection weighting matrix is W zi = σ −2 zi I N c , where I N c is an identity matrix, whose size is equal to the number of virtual measurements. The residual vectors are In the BC-DSSE context, each zero injection leads to two virtual measurements, for the real and imaginary parts of the corresponding current injection. The Jacobian rows of the ith zero injection appear as follows: where the only non-zero elements are those corresponding to the real and imaginary parts of the branch currents to and from the considered node (the sign depends on the direction assumed for the currents).
The WLS estimation step of DSSE is thus performed by means of the following N-dimensional system (where N is the number of unknowns in x) derived from (2): Possible meshes in the network can be expressed as additional constraints among the branch current state variables, using the Kirchoff voltage law, as m(x) = Mx = 0, where each couple of rows j and j + 1 refers to the Kirchoff voltage law of a mesh, expressed in real and imaginary parts, respectively. Matrix M thus contains, for each mesh, the resistances or inductances of the branches involved in the mesh itself.
Even in this case, the constraints can be included in the estimator model as virtual measurements, leading to the following: where M and W zm are the Jacobian and the weighting sub-matrix (with large weights) of the mesh constraints, respectively.

Lagrangian method
It is known that using large weights can lead to illconditioning problems in the system of normal equations to be solved. To avoid the use of very large weights, an alternative is represented by the explicit formulation of the constraints (4) in the problem. Such constrained minimization problem can be faced through the Lagrangian method [16]. According to this approach, the normal equation system (2) is extended associating at each equality constraint a Lagrange multiplier. In the case of zero injections, for a generic iteration, (2) becomes: where C is the Jacobian of the zero injection constraints, λ zi is the N c vector of the associated Lagrange multipliers and, using the same notation as in (5), G m = H T m W m H m represents the gain matrix referred only to the set of measurements and pseudo-measurements. With this approach, the number of the unknowns increases because of the multipliers, but the sparsity of system (11) is higher with respect to the case of virtual measurements.
It is worth underlining that in [15], in order to significantly reduce the numerical conditioning of the system, the use of a normalization coefficient γ = 1/ max(W ii ) for the weighting matrix is suggested. Thus, the equation system (11) has to be modified in the following form: Similarly to the zero injections, even the mesh constraints can be managed by using Lagrangian multipliers. In this case, (12) has to be rewritten as: where λ zm is the vector of Lagrange multipliers included to deal with the mesh constraints m(x).

Proposed method: state vector reduction
An efficient way to include equality constraints, wellsuited to BC-DSSE formulation, is proposed here. Each zero injection corresponds, as aforementioned, to an equivalent phasor current injection equal to zero. Since the real and imaginary parts of the branch currents are included in the state vector, it is straightforward to express the constraint on a node i as: where i is the set of branches incident to node i, while α j is a coefficient equal to +1 or −1, depending on the direction assumed for the currents in the state vector. It is easy to understand how a simple variable elimination can be performed from (14), writing one current as a function of the other incident ones. In fact, for each zero injection, a current can be expressed as a function of the other state variables and eliminated from the state. The state x can be thus written as a function of a new reduced state vectorx (of lengthÑ = N − N c ), replacing N c variables by a simple combination of the remaining state variables, as follows: where IÑ is aÑ ×Ñ identity matrix and zi is a N c × N sparse matrix of ±1 elements linking the eliminated variables x zi to the remaining ones.
Referring to the reduced state vectorx, the normal equation system (2) can be rewritten as: where the new JacobianH m is: (17) and the residual vector is computed as As for the mesh constraints, even in this case it is possible to express one of the currents involved in the mesh in terms of the remaining ones. In general, the mesh constraint is: where identifies the set of branches in the mesh, and Z j and i j are the impedance and the current of the jth branch, respectively. It is worth noting that, in a threephase framework, Z is a 3 × 3 impedance matrix, which includes also the mutual impedance terms, while i is the vector of the phase currents. From (18), it is easy to find: where k is the index of the branch whose currents will be eliminated from the state variables.
Considering this additional reduction of the state vector, it is possible to rewrite (15) in the following way: where zm is a N m ×Ñ matrix linking the N m removed current variables x zm to the remaining ones. Thus, in this case, the resulting state vector has a reduced length equal toÑ = N −N c −N m . Once the transformation matrix linking the starting state vector to the reduced one is defined, the same considerations involved in (16) and (17) hold for the execution of the estimation algorithm. The performed transformation leads to a lower sparsity of the system, reflecting the fact that each eliminated variable is expressed in terms of more remaining variables. However, since distribution grids usually have a large number of zero injections, the transformation also allows a significant reduction of the dimensions of the equation system to be solved. It is worth underlining that zi and zm are constant matrices that can be built a priori knowing the operation and topology of the network; thus, there is no need to compute them at each run of the DSSE.
It is important highlighting that such approach can be applied in an efficient way only because of the chosen state variables: a similar logic does not apply, with the same simplicity, to traditional estimators based on polar node voltages.
The proposed approach allows an efficient implementation even for the bad data detection and identification functions. In fact, the same techniques traditionally adopted in WLS estimators, based on the computation of the normalized residuals, can be conveniently implemented here, thanks to the reduced sizes of the Jacobian and gain matrices that are involved in the computation of the residual covariance matrix (see [15] for further details). Identified bad measurement data are removed, and the BC-DSSE estimation steps previously described are performed again on the reduced measurement set. It is worth noting that, in order to avoid the computation of the residual covariance matrix when bad data are not present, the bad data detection can be also implemented by using the well-known χ 2 test [15].

Test assumptions and metrics
Several tests have been performed on the unbalanced three-phase IEEE 123-bus system (Figure 1) to assess the performance of the proposed method. Data of the network can be found in [20]. Different measurement configurations have been considered in order to analyze various possible scenarios.
In order to achieve significant results from a statistical standpoint, tests have been carried out following a Monte Carlo approach. For each test, at the beginning, the true reference quantities of the network are computed by means of a load flow calculation. Then, for each trial of the Monte Carlo simulation, measurements are extracted adding random noise to the reference values, according to the uncertainty characteristics assumed for each measurement. At the end of each trial, estimation results are stored to be available for the subsequent analysis. The following assumptions are used in the tests: • number of Monte Carlo trials: N MC = 25, 000; • pseudo-measurements available for all the loads of the network and characterized by normally distributed uncertainty with maximum deviation 3σ = 50% of the nominal value; • PMU measurements characterized by uncertainty with uniform distribution and variance σ 2 equal to one third of the squared accuracy value; in particular, an accuracy equal to 0.7% and 0.7 crad (i.e. 0.7 × 10 −2 rad) is used for magnitude and phase-angle measurements, respectively, in order to simulate the accuracy limits specified in the synchrophasor standard [3].

Figure 1 IEEE 123-bus test system.
Results of the tests have been analysed to assess the performance of the proposed approach in comparison to virtual measurements and Lagrange multiplier methods. In particular, accuracy, computational properties and efficiency of the different approaches have been analysed in order to have an overall performance evaluation.

Estimation accuracy
To assess the estimation uncertainty of a given quantity, the root mean square error (RMSE) has been used. Given a generic electrical quantity y, the RMSE is defined as: whereŷ i is the estimation of y at the ith Monte Carlo trial and y true is the true value of the quantity. In this paper, as overall index for the whole network, the mean RMSE, obtained averaging the RMSEs of all the nodes or branches (depending on the considered electrical quantity), has been used.
Since the focus of the paper is on the processing of the equality constraints, a second parameter is also monitored, that is the sum of the estimated power injections in the zero injection nodes. This index allows to evaluate the capability of the method to fully satisfy the zero injection constraints. In [17], it is shown that a bad handling of the constraints in the estimator could lead to the estimation of significant power injections in the zero injection nodes, consequently affecting the overall accuracy of the estimated quantities. Equation 22 shows the definition of the indexes for both active and reactive power (with B representing the set of the zero injection nodes of the network):

Numerical properties
An important issue for the estimator design is the numerical conditioning of the equation system. In fact, as described in [15], due to ill-conditioning, small errors in the different entries of the equation system may be translated in large errors for the solution vector. As a consequence, the accuracy and the convergence properties of the algorithm can be significantly affected and numerical instabilities can appear. In all the analysed approaches, the equation systems to be solved at each iteration of the algorithm ((9), (12) and (16)) can be written in the following compact form: where G tot is a coefficient matrix, x tot is the total vector of the unknowns (note that for the Lagrangian method also the Lagrange multipliers are included in this vector) and u is a vector resulting from the measurement residuals.
To evaluate the numerical properties of the system, the condition number k of the coefficient matrix G tot is considered, that is: where · is the 2-norm of the matrix.
Other interesting properties from the computational point of view are the density and the size of the coefficient matrix. A low density (defined as the ratio between nonzero terms and total number of elements in the matrix) implies a large number of zero elements in the matrix and thus the possibility to use sparse matrix techniques for the calculations. The size, instead, is obviously associated to the number of unknowns and represents the dimension of the equation system to be solved. In the following, all these parameters will be taken into account to discuss the obtained results.

Computational efficiency
The efficiency of the estimator is a crucial factor for the real-time management and control of future distribution systems. For this reason, in the following discussion, the average execution times of the different methods (obtained averaging among the N MC Monte Carlo trials) will be compared. Moreover, since the execution times strictly depend on the number of iterations needed for the algorithm convergence, also the average number of iterations of the different approaches will be evaluated. Tests have been performed under Matlab environment and run on a 2.4-GHz quad-core processor with 8-GB RAM.

Test results
In this section, the performance of the approaches virtual measurement method (VM), Lagrangian method (LM) and state vector reduction (SVR) will be analysed and discussed.
First, a test has been carried out considering a possible realistic measurement configuration for the network. In particular, three measurement points have been supposed to be available in nodes 150 (primary station), 18 and 67. Each measurement point is composed of a voltage synchrophasor measurement on the node and of current synchrophasor measurements on all the branches converging to that node.
Results of the test show that all the analysed methods provide really similar accuracy performance. The mean RMSE of the branch active power is, as an example, 5.7 kW for all the approaches. A similar behaviour has been found also for the mean RMSEs of voltage and current. As for the estimated power injection in zero injection nodes, simulations prove that a really low total power injection can be found with a suitable setting of the algorithms, despite the large number of zero injections (119 nodes over 227 total nodes for all the three phases). In particular, both LM and SVR always give a null sum of the considered power injections. As for VM, instead, results strictly depend on the choice of the weights used for the constraints. Table 1 shows the values of P 0inj and Q 0inj achieved using different weights. It is possible to observe that increases of an order of magnitude in the weights cause corresponding decreases on the P 0inj index. Similar results have been found also for the reactive powers. In any case, it is worth noting that, since the unbalances are different in sign, the mean RMSEs basically do not change in the different scenarios. Table 2 shows the efficiency performance of the different methods. In this case, significant differences can be found depending on the approach used to deal with the zero injection constraints. In particular, it is possible to note that SVR provides the best results in terms of execution times, with an enhancement of the computational efficiency with respect to VM and LM larger than 30% and 40%, respectively. Furthermore, due to the correct modeling of the zero injection constraints, a slightly better performance can be observed for SVR even for the convergence properties (in terms of average number of iterations).
The reasons of such enhancement in the algorithm speed can be found looking at the numerical features of each approach. Table 3 reports the density properties and the size of the coefficient matrix involved in the WLS step of the estimation process. It is clear that SVR allows a large reduction of the coefficient matrix due to the elimination of all the zero injection constraints. Such reduction, despite the lower sparsity of the system, allows a faster factorization of the coefficient matrix (in this case, no sparse matrix techniques have been used for SVR) for the solution of (16). On the other hand, in LM, the introduction of the Lagrange multipliers leads to a larger size of the coefficient matrix, and the solution of the equation system, even if manageable with sparse matrix techniques, is slower.
An additional advantage guaranteed by the use of SVR concerns the condition number of the coefficient matrix. As it can be observed in Table 3, SVR gives the lowest conditioning. It is worth recalling that, in general, VM suffers possible ill-conditioning problems because of the use of very large weights to enforce the constrained measurements. As a confirmation of such an impact, the condition number obtained changing the weights to 10 8 and 10 12  9 10 10 10 11 10 12 P 0inj 1.9 · 10 −1 1.9 · 10 −2 1.9 · 10 −3 1.9 · 10 −4 Q 0inj 3.0 · 10 −1 3.0 · 10 −2 3.0 · 10 −3 3.0 · 10 −4 has been checked: in both the cases, a consequent variation of the condition number (3.74 · 10 4 and 3.02 · 10 8 , respectively) has been found. Several tests have been performed also to verify the proper operation of the aforementioned bad data detection and identification function when using the proposed approach. For instance, tests have been performed by adding intentional errors (of 5%, 10% and 20%) to the voltage magnitude measurement at node 18 (on the first phase). In such a scenario, all the analysed methods allow the detection of the bad data, through the χ 2 test, and the proper identification of the erroneous voltage, by means of the largest normalized residual technique. It is important to note that, since the presence of the bad data implies the computation of the residual covariances and the need to repeat the estimation process, the aforementioned improvements brought by the SVR method on the execution times further increase. It is also worth noting that, in some cases, depending on measurement configuration, when considering the bad data on other measurements (for example, on PMU currents), the bad data identification function could be unable to properly identify the erroneous measurement. However, as known from the literature, this is a general issue that can occur due to the low redundancy of the measurements, and it does not depend on the particular method used to handle the equality constraints. In fact, several tests (not reported here for the sake of brevity) have been performed changing the corrupted measurement, and as expected, test results prove that all the approaches exhibit exactly the same behaviour, with the same identification results, in all the different scenarios. Since possible identification problems are generally related to the measurement system deployed in the distribution grid, a deeper analysis on this issue is out of the scope of this paper.

Impact of the measurement configuration
Further tests have been performed to assess the performance of the methods with different measurement types and configurations. First of all, the general validity of the previous considerations has been tested using different measurement devices. To this purpose, voltage and current phasor measurements have been replaced by voltage magnitude and active and reactive power measurements (with accuracy equal to 1% and 3%, respectively). It is worth noting that in this case, since synchronized measurements are not available, the estimator model has to be suitably adapted to take into account the absence of an absolute phase angle reference (see [12] for more details). Tables 4 and 5 show the results about the efficiency performance and the numerical properties of the different methods in this case with traditional measurements. It is possible to observe that all the previously made considerations still hold, with SVR providing the best performance in terms of execution time and numerical conditioning. As for the accuracy performance, even in this case, all the methods provide really similar results.
As for the measurement configuration, the attention has been mainly focused on the impact of additional voltage measurements, since they can significantly affect numerical properties and efficiency of the branch current-based formulation of DSSE. In fact, voltage measurements lead to non-zero terms in the Jacobian H m corresponding to all the derivatives with respect to the branch currents included in the path between the bus used as reference in the state vector and the measured node (for details, see [12]). This, in turn, causes a lower sparsity of the coefficient matrix in (23), thus affecting the efficiency of the equation system solution.
As an example, Tables 6 and 7 show the results about the numerical properties and the computational efficiency when two voltage measurements (in nodes 86 and 105) are added with respect to the configuration used in the previous tests.
As expected, the presence of additional voltage measurements affects the coefficient matrix leading to higher densities: this is the main reason for the increased execution times shown in Table 7 (with respect to the results in Table 2). Moreover, it is worth highlighting that the increased density brings different impacts on the different  methods. In fact, since the solution of the equation system in SVR is managed without using sparse matrix techniques, the impact on this approach is smaller with respect to the other methods. This is confirmed by the enhancements obtained on the execution times that, in this scenario, rise up to almost 37% and 50% with respect to VM and LM, respectively. To emphasize and check such result, a further test has been carried out adding four supplementary voltage measurements (in nodes 25, 42, 48 and 91). Table 8 shows the results concerning coefficient matrix density and execution times. It is possible to observe that further increases in the density of the coefficient matrix, with consequent further improvement of the computational efficiency of SVR, are achieved. In this case, the time saved through SVR is larger than 45% with respect to VM and 60% with respect to LM.
As for the accuracy and the numerical conditioning of the analysed methods, considerations similar to those made for the first test can be derived also in these cases: all the methods provide very similar accuracy results, and SVR shows the best conditioning properties.

Impact of the size of the network
One of the main issues involved in the handling of distribution systems is the large size of these networks. This aspect is particularly critical from the standpoint of the execution times, since it implies a significant increase of the size of the equation system to be solved. In such a situation, possible drawbacks can arise for the SVR method due to the fill-ins resulting from the elimination of the state variables.
To assess the impact of such issue, additional tests have been performed with networks having a larger number of nodes. In particular, in order to simulate different sizes   Figure 1). In order to obtain different loading conditions in the feeders, power consumptions of the loads have been modified adding random variations. Performed tests show an obvious decrease for the density of the coefficient matrices in all the tested methods. For instance, in the SVR approach, the coefficient matrix density drops to 4.42% and to 2.22% when five and ten feeders are considered, respectively. In such a situation, even considering the huge size of the coefficient matrix, the use of sparse matrix techniques is a forced choice also for the SVR approach. Even in this case, the SVR method maintains the already cited advantages in terms of execution times. As an example, Table 9 shows the obtained results, for iteration numbers and execution times, when ten feeders are considered in the network (the resulting three-phase grid is, in this case, composed of more than 2,000 nodes). It is possible to observe that, because of the large size of the equation system to be solved, the execution times are now significantly higher than the previous cases. However, despite of the increased size of the network, the SVR method still allows a significant enhancement of the computational performance, leading to an improvement of about 22% and 37% with respect to the VM and the LM approaches, respectively.

Impact of the mesh constraints
In this section, the performance of the proposed method are tested for the case of weakly meshed networks,  referring to the original 123-bus grid. To this purpose, the presence of a branch between nodes 151 and 300 and between nodes 54 and 94 of the benchmark network has been supposed in order to create two meshes. Simulations have been carried out considering the base monitoring configuration composed of the three measurement points in nodes 150, 18 and 67. Table 10 shows the numerical properties of the coefficient matrix for the different methods.
It is possible to observe that the density of the coefficient matrix in VM and SVR is significantly affected by the presence of the meshes. In fact, as clear from (18), the mesh constraint involves all the three-phase currents of the branches belonging to the mesh. Thus, the matrix multiplications needed to create the coefficient matrix G tot and involving the Jacobian matrix M for VM (see Equation 10) and the transformation matrix zm for SVR (see Equations 20,16 and 17) lead to a significant increase of the non-zero elements. At the same time, the presence of the additional branch currents (due to the meshes) lead to an increase of the dimension of the equation system for VM, while no change appears for SVR since such currents are expressed in terms of the reduced state vector. In the case of LM, instead, the explicit expression of the mesh constraints, together with the growth of the dimensions of the equation system (due to both the additional branch currents and the additional constraints) allows to keep the coefficient matrix very sparse.
All these aspects bring direct effects on the efficiency of the different methods. Table 11 reports the obtained results for average iteration number and execution times. It is possible to observe that in this scenario, the VM approach is negatively affected by the presence of the constraints and gives the worst results in terms of execution time. Instead, the proposed method still shows significant  benefits on the computational performance with improvements on the execution times around 45%. The same considerations already made for the tests with the radial version of the network hold for the accuracy performance and the conditioning properties.

Conclusions
In this paper, an efficient way to include the equality constraints in a branch current-based state estimator is presented. The method exploits the use of rectangular branch currents as state variables of the system to perform a simple elimination of one of the currents involved in a constraint, expressing it as a linear function of the remaining ones. The method not only is particularly efficient in the management of zero injections but also allows the treatment of mesh constraints. Performed tests prove the goodness of the proposed technique and in particular its capability to significantly improve the computational efficiency of the estimator (with respect to other traditionally used methods). Moreover, full fulfilment of the constraints is guaranteed, and additional benefits can be achieved for the numerical conditioning of the system.