Multi Objective for PMU Placement in Compressed Distribution Network Considering Cost and Accuracy of State Estimation

A phasor measurement unit (PMU) can provide phasor measurements to the distribution network to improve observability. Based on pre-configuration and existing measurements, a network compression method is proposed to reduce PMU candidate locations. Taking the minimum number of PMUs and the lowest state estimation error as the objective functions and taking full observability of distribution network as the constraint, a multi objective model of optimal PMU placement (OPP) is proposed. A hybrid state estimator based on supervisory control and data acquisition (SCADA) and PMU measurements is proposed. To reduce the number of PMUs required for full observability, SCADA measurement data are also considered into the constraint by update and equivalent. In addition, a non-dominated sorting genetic algorithm-II (NSGA-II) is applied to solve the model to get the Pareto set. Finally, the optimal solution is selected from the Pareto set by the technique for order preference by similarity to ideal solution (TOPSIS). The effectiveness of the proposed method is verified by IEEE standard bus systems.


Introduction
High penetration of distributed generator and controllable load has changed the fundamental way of the operation in the distribution network from a passive to active network [1][2][3].The distribution network is developing in an intelligent direction [4].The uncertainty of distributed generator and load brings great challenges to the stable operation of the distribution network, such as bidirectional power flow and increased voltage fluctuation.Furthermore, the distribution network topology changes dynamically to reduce power loss and voltage deviation [5].As uncertainty increases, real-time measurement devices need to be added to the distribution network.Improving the observability of the distribution network is an important way to realize the intellectualization of the distribution network.Monitoring network state is traditionally done through the supervisory control and data acquisition (SCADA) system whose measurements are usually unsynchronized [6].Among those newly-developed devices that aim at realizing the online system monitoring, phasor measurement unit (PMU) opens a new avenue for the distribution system.PMU is a real-time measurement device based on global positioning system (GPS), and it can receive the GPS timing signal to realize the synchronous acquisition of bus data.Compared with other measurement devices, PMU can deliver highly synchronized direct phasor measurements because it has high GPS synchronization accuracy.Therefore, the PMU can not only be applied to real-time state estimation of the distribution system and optimize its operating state, but also can quickly deal with problems such as anomaly detection [7], overcurrent protection [8], steady-state model synthesis [9], network topology detection [10], reliability analysis [11] and so on.Reference [7] proposes a layered architecture for monitoring networks based on the PMU measurements and establishes a set of analyses and sensor primitives to detect anomalous behavior in the control perimeter.Reference [8] proposes a scheme to enhance the traditional overcurrent protection scheme by using PMU functions embedded in solid state transformer and feeder protection relay.The proposed solution not only can enable the relay to operate correctly at its original range but also can maintain coordination with each other.In Reference [9], a new method of synthesizing the steady-state model of the unbalanced distribution network is proposed, which uses real-time measurements provided by PMU and can reduce the complexity of the network model caused by the increase of the power supply in the distribution network.In Reference [10], a new method for detecting topology changes in the distribution network is proposed, which is called the time series signature verification method for topology detection.This method analyzes the measurements provided by PMU installed on the distribution feeder.Reference [11] analyzes the reliability of the distribution system in the presence of the adjusted download resources with the PMU measurements from the same distribution feeder, which considers the practical regulation-eligible load types, the practical distribution-level protection devices, the impact of delay and so on.Therefore, PMU is considered as one of the most important measurement devices, and utilization of PMU in the distribution network is of great interest.
PMU can acquire the voltage phasor of the bus and the associated branch current phasor at high frequency [12].When the voltage phasor of a bus can be measured directly or indirectly, the bus is called an observable bus.When all the buses of the distribution network are observable, the network is called an observable network.If PMU is installed on each bus, the network is fully observable.However, the installation cost is high, and it is difficult to achieve the full deployment of PMU in the short term.How to achieve optimal PMU placement (OPP) in a distribution network encounters technical challenges.
The purpose of OPP is to find out the most representative location and place a PMU on it.Around the theme of improving the level of observability, the main goals of OPP include minimizing the number of PMUs required, maximizing measurement redundancy and handing contingency constraint.The main goal is to minimize the number of PMUs required.In Reference [13], an OPP method for complete observability of a system under normal operating conditions is proposed.Tabu search algorithm is used to determine the minimum number of PMUs required to make the system observable.A new PMU layout formula is developed to ensure the full observability of the network in Reference [14].It takes the network observability under all possible arrangements of line connections at complex buses into account.Another goal is to maximize measurement redundancy.Measurement redundancy represents the number of times the bus is observed directly or indirectly.The solution with the largest network measurement redundancy is chosen as the best choice in Reference [15].Network measurement redundancy is the sum of the measurement redundancy of all buses.In Reference [16], a multi-objective probability model is proposed for the problem of OPP.The proposed model simultaneously optimizes two objective functions: minimizing the number of PMUs and maximizing the expected value of network measurement redundancy.Various accidents may occur in the network.Therefore, it is important to analyze the observability in these conditions which include line outage and PMU failure.In Reference [17], a new OPP method which considers emergency constraints is provided.The proposed method is based on n-k redundancy criterion using robust optimization.This security standard ensures network observability in any accidental state, including up to k loss of PMUs.In Reference [18], the OPP problem is performed in order to achieve a fully observable network under normal and emergency conditions considering network expansion.To this end, the strain index is introduced through the modeling process of n-1 contingency states.
Unlike the transmission grid, the measurement data of the distribution network is insufficient, which renders most distribution networks unobservable [19].Existing real-time measurements of the distribution network are mainly provided by the SCADA system [20], and the remote terminal unit (RTU) is the main measurement device.On the one hand, the sampling period of SCADA is much longer than that of PMU [21].On the other hand, the SCADA system cannot provide the phase angle data directly.The lack of a sufficient number of measurements is usually solved by introducing pseudo-measurements, which are mainly represented by historical measurements or forecast measurements on generators production and power consumption [22,23].
The main purpose of improving network measurement redundancy is to improve the accuracy of state estimation.The distribution network state estimation technology is expected to be able to monitor the operating status of the system accurately and provide reliable data for other advanced management software of the active distribution network, such as transient stability [24], voltage security visualization [25], fault location identification [26], dynamic control strategies [27] and so on.Therefore, taking state estimation accuracy instead of network measurement redundancy into the OPP problem can be more helpful to improve the observability and controllability of the distribution network.Distribution system state estimation methods proposed are mainly based on weighted least squares (WLSs) algorithms [28][29][30], and considerable effort has been devoted to reducing computational requirements.The main differences between proposals are basically the choice of state variables, the simplification of accelerated estimation and the technique of combining heterogeneous measurements.The two main kinds of state variables are bus voltages and branch currents.The weight associated with measurement is proportional to the accuracy of the measurement.Reference [28] discusses the correlation and importance of the weighted least squares estimation approach and presents the impact of the input data correlation on state estimation results.The implementation details of the constrained weighted least squares state estimation method are proposed in Reference [29], which can implement the equality constraint formula of the WLS method.Reference [30] uses real-time measurements and all available pseudo-measurements, the model is extended to include the slack bus voltage during the estimation process, and the possibility to treat radial and weakly mesh topologies is also present.
In order to reduce the number of PMUs required for the full observability, zero injection bus (ZIB) should be considered into the OPP model.If a bus is not connected with any load, and the current injection of the bus and the power measurements of the bus are both zero, this kind of bus can be called a ZIB [6].The observability of the distribution network can be improved by the presence of ZIB by Kirchhoff laws [31].The existing methods of taking ZIBs into consideration are mainly divided into three types: linear inequality constraint, nonlinear modification and topology transformation.Reference [17] incorporates ZIBs to the OPP model by adding a linear inequality constraint.The number of added inequality constraints is equal to the number of ZIBs.Although this kind of method can guarantee that the model is nonlinear, this kind of method is difficult to apply when there are other measurements or multiple ZIBs are connected in the system.Topology transformation refers to merging the ZIB with one of the connected buses.Nonlinear modification means that the ZIB is applied to the OPP problem by modifying existing constraints.This method is easy to implement and has good adaptability, but it also makes the model nonlinear.In Reference [32], selection rules for topology transformation are proposed.ZIBs are merged with one of their connected buses.Based on the three rules, the proposed method determines the best bus to merge with the ZIB to get the minimum number of PMUs required for the full observability.However, using topology transformation to deal with ZIBs has two limitations: to identify the locations of PMUs and the importance of selecting the merged bus.Although the number of PMUs required for the full observability is reduced, the optimization effect is affected by the selection of the bus which is merged with ZIB.
Due to a large number of buses in the distribution network, only using PMU measurements to realize the network observable fully is not economic.Therefore, conventional measurements should also be applied to the OPP problem.Conventional measurements of the distribution network are mainly provided by the SCADA system.In Reference [33], conventional measurements are modeled, and a recipe is introduced when multiple branches are adjacent to the same bus.The changes in the flow measurements after branch interruption are also considered in the proposed OPP method to obtain a more practical solution for the distribution network.Only power flow measurements are used in this OPP model, and other useful measurements are not utilized.Reference [34] presents a new integration model which considers the impact of current measurement and ZIB on the OPP model.The OPP model considering current measurements and ZIB is based on the logical sum of the observability of the bus, which can take advantage of the maximum capabilities of current measurement and ZIB.In addition, this method can also be applied to other systems with multiple topologies.
Compared with the transmission network, the number of buses in the distribution network is large, so it is difficult to solve the model and verify the optimal solution.In order to search the optimal solution efficiently, some OPP methods are proposed to reduce the PMU candidate locations and improve the optimization effect.In Reference [35], considering the radiation of the distribution network, an OPP method which decomposes the network into two sub-networks is proposed.This method decomposes the whole network into two sub-networks: the primary sub-network and secondary sub-network.The OPP process of each sub-network is performed separately, so the number of PMU candidate locations is greatly reduced.Although the proposed OPP method saves calculation time, the number of PMUs required for the full observability increases due to the segmentation operation.
Numerical algorithms [36,37] and heuristic algorithms [38][39][40][41][42] are the main optimization techniques for solving the OPP model.The numerical algorithm is mainly applicable to solve the problem whose model is simple and linear.Reference [36] proposes an iterative weighted least squares algorithm to solve the OPP model for the full observability, which has actual placement variables.The proposed algorithm is implemented by using the available routines in the existing state estimator to avoid the use of other mathematical optimization technologies.Reference [37] proposes a mathematical OPP framework which considers cost and benefits at the same time.The enforcement in reliability is calculated by using the PMU measurement infrastructure and converted to a monetary value as the development benefit.The benefit is reduced by the cost value to get the optimal solution.The heuristic algorithm is applicable to solve the nonlinear, high-dimensional and multi-objective model.In Reference [38], a greedy algorithm is applied to the optimization.The proposed greedy algorithm has good performance with low computational complexity and robustness to solve the uncertainty of PMU placement budget.This performance makes the greedy algorithm attractive in typical scenarios of phased PMU placement because the performance is robust to changes in the PMU budget.Reference [39] develops an effective exponential binary particle swarm optimization algorithm to obtain multiple solutions to OPP.An exponential inertia weighting factor is introduced by the proposed algorithm to improve the group's searchability.In order to combine the previous positions of the particles, two innovative equations are presented to update the particle position.In addition, two useful filtering techniques are applied to converge reliably.Reference [40] proposes a biogeography-based optimization algorithm to solve the OPP problem.This proposed algorithm is suitable to solve the multi-objective optimal PMU placement problem for the full observability.In Reference [41], a backtracking search algorithm is applied to the OPP problem.Selection operation, variation operation and crossover operation of backtracking search algorithm are different from those of other algorithms.The backtracking search algorithm can create more effective individuals in each generation.It controls the magnitude of the search direction in a balanced way that is required for global and local searches.Reference [42] proposes a practical multi-stage simulated annealing algorithm for the optimization of PMU placement.The proposed OPP method can also solve the problem of PMU channel limits.In order to control uphill and downhill motion, the multi-stage simulated annealing algorithm applies multiple stages in the optimization process, which can encourage to converge to the OPP more quickly than classical simulated annealing.
In this paper, a new OPP method for a distribution network is proposed, which considers cost and state estimation accuracy.The main contributions are: (1) Based on pre-configuration and existing measurements, a network compression method is proposed, which can reduce the candidate locations of the PMUs and improve the network measurement redundancy.Additionally, a hybrid state estimator is applied, which can make full use of supervisory control and data acquisition (SCADA) and PMU measurements to improve state estimation accuracy.
(2) A new OPP model is built which can minimize the number of PMUs required for the full observability and state estimation errors of the distribution network simultaneously.Based on this, it can maximize network measurement redundancy.
(3) To reduce the number of PMUs required for the full observability, a new way of taking SCADA data into the constraints is also developed by updating some constraints and equating some buses to "lesser ZIB."To solve the OPP model, a non-dominated sorting genetic algorithm-II (NSGA-II) is applied to get the Pareto set, and the compromise solution is selected by the technique for order preference by similarity to ideal solution (TOPSIS) which can comprehensively minimize cost and state estimation errors.The optimal solution which has the minimum number of PMUs to ensure the network observable can also be provided by the Pareto set.
The remainder of this paper is organized as follow.Section 2 proposes the network compression method.Section 3 presents the hybrid state estimator which considers PMU measurements.Section 4 proposes the multi-objective model of OPP.Section 5 provides the solution methodology of the model based on NSGA-II and TOPSIS.Section 6 presents the simulation results and Section 7 concludes the paper.

Network Compression
Network compression includes two parts.On the one hand, some observable buses are merged with their neighbors by pre-configuring a part of PMUs; on the other hand, based on existing ZIB measurements, some buses which can become an observable bus are merged with their neighbors.Network compression can reduce the PMU candidate locations, increase network measurement redundancy and improve the optimization effect.

The Observability of Distribution Network
For the distribution network, when the voltage phasor of a bus can be measured directly or indirectly, the bus is identified as observable.PMU can acquire the voltage phasor at the installed bus and the current phasor of all the branches connected to the installed bus at high frequency.Based on the characteristics of PMU, the rules of observability are as follows: (1) If a bus is placed with a PMU, this bus and its neighbors are identified as observable buses (neighbors refer to buses which are connected to the bus).
(2) If all neighbors of ZIB are observable buses, the voltage phasor of this ZIB can be calculated by Kirchhoff laws (KCLs), so this ZIB is also identified as an observable bus.
ZIB has no load or generator.It acts as a transshipment bus in distribution networks.The sum of currents flowing into a ZIB is zero.
(3) If a ZIB is an observable bus and the observability of only one of its neighbors is unknown, the voltage phasor of this neighbor can be calculated by KCLs, so this neighbor becomes an observable bus.

Network Compression Based on Pre-Configuration and Existing Measurements
When the distribution network has n buses, the PMU has n candidate locations.Based on the rules of observability, the bus which needs to be placed with PMU and the buses which do not need to be placed with PMU can be determined.Through pre-configuration, the bus which does not need to be placed with PMU is merged with its neighbor to reduce PMU candidate locations.The specific methods are as follows: (1) Place a PMU on the neighbor of the leaf bus, and merge the leaf bus with its neighbor into a new bus.
The leaf bus is a bus which is only connected to one bus.The distribution network is radial, so it always has many leaf buses.When a PMU is installed on a leaf bus, the number of new observable buses is only two: leaf bus and its neighbor; when a PMU is placed on the bus which is connected to a leaf bus, more than two buses become observable buses.Therefore, PMU is not placed on the leaf bus, PMU is placed on the bus which is connected to the leaf bus instead.Therefore, the leaf bus becomes an observable bus, which does not need to be considered in the process of OPP.The leaf bus can be merged with its neighbor.Network compression based on pre-configuration can reduce the number of buses in the network and improve network measurement redundancy when there are many leaf buses in the distribution network.
(2) Merge the ZIB with one of its neighbors while only two buses are connected to the ZIB.
For the ZIB which has only two neighbors: (a) if a PMU is placed on any one of the neighbors, both the ZIB and another neighbor become observable; (b) if all two neighbors of the ZIB are observable buses, the voltage phasor of the ZIB can be calculated by KCLs, so this ZIB is identified as an observable bus.Therefore, if the ZIB has only two neighbors, this ZIB does not need to be placed with a PMU and can be merged with one of the neighbors.
In summary, the network compression method provides a compressed network for the OPP model based on the pre-configuration and existing measurements.

The Basic Model of State Estimation
The state estimator proposed in this paper is based on WLS.Assuming that the number of state variables is N and the number of measurements is M, the measurements can be expressed as Equation ( 1): where T is the M-dimensional vector of measurement errors, whose entry v i represents the error of measurement i; h is the vector of functions, which represents the relationship between z and x.According to the measurement vector z and measurement function h, the objective function that minimizes the sum of squares of estimation errors can be expressed as Equation ( 2): where J(x) is the sum of squares of estimation errors; W is the weight matrix of measurements, diag W = [w 1 , w 2 , • • • , w M ] T ; w i represents the weight of measurement i, which can be described as Equation ( 3): where R is the covariance matrix of the measurement error; σ i is the variance of the measurement i.
To minimize the sum of squares of estimation errors, J(x) is guided by x, which can be expressed as Equation (4): According to Equation (2), Equation ( 4) can be modified to Equation ( 5): where H(x) is named as the Jacobian matrix, which can also be defined as Equation ( 6): As h is the vector of nonlinear functions, the Taylor series is applied to solve Equation ( 5), and the high terms with orders greater than 1 are ignored to get a linear Equation (7): where x 0 represents the initial state vector.According to Equations ( 5) and ( 7), the iterative equations can be expressed as Equation ( 8): where x (k+1) and x (k) are the state variable vectors obtained by the k − th and (k + 1) − th iterations respectively.Newton method is applied to solve the iterative Equations ( 8).The estimated state vector x can be obtained when the convergence criterion is satisfied.The convergence criterion can be defined as Equation (9): where ε represents the value of convergence.

Estimator with PMU Measurements
The state estimator based on SCADA measurements is nonlinear, while the state estimator based on PMU measurements can be linear.In this paper, a hybrid state estimator considering PMU measurements and SCADA measurements is proposed.
At the sampling time of SCADA, the first state estimator is carried out with the SCADA measurements and PMU measurements at the nearest moment, and the result of state estimation xs is saved.This state estimator is nonlinear.At this moment, the measurements include the voltage amplitude and phase angle of the bus, the amplitude and phase angle of branch current, the active power and reactive power of branch, the active power and reactive power injected at the bus.In this paper, the real and imaginary parts of slack bus voltage and branch current are chosen as state variables [30].The measurement function and Jacobian matrix are obtained by the relationship between the measurements and state variables.PMU collects the amplitudes and phase angles of bus voltage and branch current, while the state variables of the estimator are the real and imaginary parts of bus voltage and branch current, so PMU measurements need to be equivalently transformed.Coordinate transformation involves the transmission of measurement errors.Taking the branch current as an example, the covariance of branch current in terms of real and imaginary parts can be calculated as Equation ( 10): where w and k are the bus number; σ 2 are the measurement error covariances of amplitude and phase angle of branch current.Similarly, the measurement errors of the real and imaginary parts of the slack bus voltage can be obtained.Then, at the sampling time of PMU, the second linear state estimator is carried out with the last result of state estimation xS and PMU measurements at the current moment.The covariance matrix of the estimated state error can be calculated as Equation (11): where R S represents the covariance matrix of the estimated state error of xS .
Then the weight matrix of measurements of this state estimator can be expressed as Equation ( 12): where W P is the weight matrix of measurements of this state estimator; R I PMU and R U PMU are the covariances of branch current and bus voltage measurement errors provided by PMU at the current moment.The estimated state provided by the second linear state estimator can be defined as Equation ( 13): where xP is the estimated state provided by the second state estimator; H P is the Jacobian matrix of this state estimator, its elements are constant which can be obtained by Equation (6).
In this paper, the state estimation error is provided by the error in voltage amplitude estimation and error in voltage phase angle estimation.
The error in voltage amplitude estimation can be expressed as Equation ( 14): where n is the number of bus in the network; e amp is the mean error in voltage amplitude estimation at each bus, which represents the error in voltage amplitude estimation; U i,true is the true value of bus voltage amplitude provided by the power flow computation; U i,se is the estimated value of bus voltage amplitude provided by the proposed hybrid state estimator.The error in voltage phase angle estimation can be expressed as Equation ( 15): where e ang is the mean error in voltage phase angle estimation at each bus, which represents the error in voltage phase angle estimation; ϕ i,true is the true value of bus voltage phase angle provided by the power flow computation; ϕ i,se is the estimated value of bus voltage phase angle provided by the proposed hybrid state estimator.

Multi Objective Model of Optimal PMU Placement
Considering the cost and accuracy of state estimation, a new multi-objective model of optimal PMU placement is presented.The objective functions are to minimize the number of PMUs and minimize the error of state estimation simultaneously, and the constraint is to ensure the distribution network is observable.On the basis of minimizing the number of PMUs required and the state estimation error, the solution with the largest measurement redundancy is obtained by modifying the objective function, so as to deal with the problem of data loss caused by communication problems.
In addition, SCADA measurements are also applied to reduce the number of PMUs required for the full observability.

The Basic Model of OPP Considering Cost and State Estimation Accuracy
The first goal is to minimize the cost, so the objective function f 1 (X) is the number of PMUs, which can be expressed as Equation ( 16): where N S is the number of buses in the compressed network; i is the bus number; X = is the N S -dimensional vector of placement; X i is a 0-1 variable, which can be expressed as Equation ( 17): The second optimization objective is to minimize the state estimation error, that is, to maximize the state estimation accuracy.The objective function f 2 (X) is the sum of the error in voltage amplitude estimation and error in voltage phase angle estimation, which can be expressed as Equation (18): The constraint is that the distribution network is fully observable, which can be expressed as Equation (19): where M is the measurement redundancy vector, which presents measurement redundancy of each bus, and its entry m i is the number of times that bus i is observed directly or indirectly T is the minimum measurement redundancy vector, b i = 1 ensures each bus observable; A is the connectivity matrix, which can be defined as Equation (20): 1 if bus i is connected to bus j The more buses the network has, the more solutions with the same number of PMUs required for the full observability exist.Therefore, the network measurement redundancy becomes an important index to select the optimal solution [43].In order to ensure that the optimization effect of searching the minimum number of PMUs does not be affected, the first objective function f 1 (X) can be modified to Equation (21):

Updated Constraint Considering ZIB Measurement
Assuming that a ZIB is connected to q buses, and any q of these buses are observable, the voltage phasor of us (q + 1) can be calculated by KCLs.In other words, bus (q + 1) is also observable.If a ZIB is connected to two buses, this ZIB can be merged with one of its neighbors.When the network has some ZIBs which have more than two neighbors individually, the Equation ( 19) can be updated to Equation (22): where is the ZIB measurement vector, which represents whether ZIB and its neighbor can be observed by the characteristic of ZIB, its entry o i can be defined as Equation (23): where g i can be calculated as Equation ( 24): where bus p is a ZIB, from bus (p + 1) to bus (p + q) are its neighbors.

Updated Constraint Considering SCADA Measurement
The installation cost can be high if only PMU measurements are applied to meet the constraints.In the process of OPP, the existing SCADA measurements are applied to meet the full observability.
The full observability indicates that the voltage phasors of all buses in the network are known.In this paper, power flow measurements and injection measurements are applied to the OPP model by updating the constraints and equivalent.
The relationship between power flow and bus voltage phasor can be expressed as Equation ( 25): where S ij is the power flow of branch i − j; • U i represent the voltage phasor of bus i; * U i and * U j represent the conjugations of the voltage phasor of bus i and bus j respectively; * y ij represents the conjugation of branch admittance.
According to Equation ( 25), if the power flow and the voltage phasor of the bus at one end of the branch are known, the voltage phasor of the bus at another end can be calculated.
If the power flow of branch i − j is measured, the constraints corresponding to bus i and bus j can be modified to Equation ( 26): where A i and A j are respectively the i − th and j − th row of A; k i and k j are 0-1 variables, which can be calculated as Equation ( 27): If k i = 1, the bus i is observable.If k i = 0, the bus i is unobservable.
The relationship between the power injection and the bus voltage phasor can be expressed as Equation ( 28): where S i represents the power injection at bus i; * I i represents the conjugation of the current phasor of the bus i; j ∈ i represents bus j is connected to bus i in the network; * Y ij represents the conjugation of the admittance matrix element.
According to Equation ( 28), if the power injection at the bus i is known, and the number of observable buses (bus whose voltage phasor is known) among bus i and its k neighbors is k, the voltage phasor of bus (k + 1) can be calculated by KCLs.The observation property of the bus whose injection is known is the same as that of ZIB.
Therefore, if the power injection at the bus is known, this bus is equivalent to ZIB, which is named as "lesser ZIB" in this paper.If the "lesser ZIB" has two neighbors, the "lesser ZIB" is merged with one of its neighbors and the network is compressed again.If the "lesser ZIB" has more than two neighbors, the modeling process is shown in Equations ( 22)- (24).
The reason why this kind of bus is named as "lesser ZIB" is that ZIB measurement is perfect while the power injection measurement has the errors.Although the bus with known power injection has some similar characteristics of ZIB, the weights should be different when these measurements are applied to state estimation.

The Final OPP Model
To sum up, the multi-objective model of OPP considering cost and state estimation accuracy can be expressed as Equation ( 29): (29)

Solution Methodology Based on NSGA-II and TOPSIS
The proposed OPP model is nonlinear because SCADA measurements and ZIBs are introduced into the problem.The multi-objective model is usually transformed into a single-objective model by unifying the dimensions and setting the corresponding weights [44].However, the selection of weight requires prior knowledge which has strong subjectivity.Therefore, the non-dominated sorting genetic algorithm-II (NSGA-II) is employed to perform the optimization and get the Pareto optimal solution set in this paper.
To minimize the number of PMUs required and state estimation error comprehensively, the compromise solution is selected by TOPSIS [45].If the installation cost is limited, the solution with the minimum number of PMUs in the Pareto set is selected as the optimal solution.

Initialization Based on Network Compression
The network is compressed by pre-configuration and existing measurements.The buses in the compressed network should be renumbered according to the ascending order of the original number.
The individuals are encoded in binary, and the length of the string is N S .If bus i is placed with PMU, the corresponding gene is set as 1.If bus i is not placed with PMU, the corresponding gene is set as 0. During the process of generating the initial population, the corresponding genes of the buses pre-configured with PMUs are uniformly set as 1, while the corresponding genes of other buses are randomly set as 1 or 0.

Fitness Calculation
Combining two objective functions with the constraint in the form of the penalty function, the fitness functions can be expressed as Equation (30): where F 1 (X) and F 2 (X) are fitness functions of f 1 (X) and f 2 (X) respectively; C is the penalty factor, whose value is much larger than 1 to reduce the fitness greatly if the solution does not achieve full observability.

Genetic Operations Based on Pre-Configuration
In this paper, the probabilities of crossover operation and mutation operation are no longer set to constants.Instead the probabilities of crossover and mutation change with the fitness of the individual.The specific process of calculating the probabilities of crossover and mutation is shown in Reference [46], which is no longer described here.
The corresponding genes of the buses pre-configured with PMUs should not change.Therefore, after the crossover operation and mutation operation, the corresponding genes of the buses pre-configured with PMUs are re-assigned to 1.

Fast Nondominated Sorting
If individual X i is not dominated by individual X j , the relationship between them can be defined as Equation (31): An individual can be called as a non-dominant individual if it is not dominated by any other individuals in the population.According to the fast nondominated sorting approach, the group is stratified, and each individual gets the corresponding rank i rank .

Crowding-Distance Calculation
According to the value of F 1 (X), the individuals of each rank are rearranged in ascending order, the crowding distance can be calculated as Equation (32): where T i represent the crowding distance of individual i.
According to the value of F 2 (X), the individuals of each rank are rearranged in ascending order again, and the crowding distance can be updated to Equation (33): The crowding distances of individuals at two ends of each rank are set as a large value, so individuals at two ends have an advantage in selection.
In the process of optimization, the individual with small rank is superior to the individual with large rank.For the same rank, the individual with a large value of crowding distance is superior to the individual with a small value of crowding distance.

Optimization Process
The steps of applying the above NSGA-II to solve the proposed OPP model are as follows: (1) Input the information of distribution network: the connectivity matrix A, locations of ZIBs, SCADA measurements and so on; (2) Compress the network and rearrange bus number; (3) Create the initial population which has Q individuals and choose this population as the parent population; (4) Set the total number of iterations MaxGen and set iteration number as Gen = 1; (5) Perform crossover operation and mutation operation to generate the offspring population, and merge the offspring population with the parent population; (6) Calculate the rank and crowding distance of each individual; (7) Select Q individuals according to the rank and crowding distance to generate the parent population, and set Gen = Gen + 1; (8) Whether the maximum number of iterations is reached?Yes, end; No, go to step (5).

Selection Based on TOPSIS
The technique for order preference by similarity to ideal solution (TOPSIS) is a kind of comprehensive analysis method according to the distance between the evaluation object and the ideal solution.In this paper, TOPSIS is applied to select the compromise solution from the Pareto set.The process of TOPSIS algorithm is as follows: At first, the values of f 1 (X) and f 2 (X) are applied to construct a decision matrix, which can be defined as Equation (34): where T is the decision matrix, and t ij represents the decision indicator obtained by the value of the objective function j of the solution i; Q is the number of solutions in the Pareto set.By normalizing T, the standard matrix P is obtained, which can be expressed as Equation (35): Then, the entropy values of the two objective functions can be calculated as Equation ( 36): where e j represents the entropy value of the objective function f j (X).
According to the e j , the entropy weights of two objective functions can be calculated as Equation (37): (1 where α j represents the entropy weight of the objective function f j (X).
Considering the subjective weight of the planner, the weights of two objective functions can be calculated as Equation ( 38): where w j represents the final weight of objective function f j (X); λ j is the subjective weight of the planner, the larger λ j is, the more important f j (X) is.According to w j and P, the weighted decision matrix Z is constructed, which can be defined as Equation (39): Next, the extreme values of the weighted decision matrix can be calculated as Equation (40): where o + j and o − j represent the maximum value and minimum value of the j − th column of the weighted decision matrix Z respectively.
After that, the positive ideal solution and the negative ideal solution can be calculated as Equation (41): where Y + and Y − represent the positive ideal solution and the negative ideal solution respectively.For each solution in the set, the distances from the positive ideal solution and the negative ideal solution can be calculated as Equation (42): where d + i represents the distance between the positive ideal solution and solution i; d − i represents the distance between the negative ideal solution and solution i.
Then, the closeness coefficient of each solution in the Pareto set can be calculated as Equation ( 43): where D + i represents the closeness coefficient of solution i.At last, the solution with the largest value of the closeness coefficient D + i is selected as the compromise solution.
The flow chart of the multi-objective OPP method in the compressed distribution network considering cost and state estimation accuracy is shown in Figure 1.Calculate errors in voltage estimation by equations ( 14) and ( 15)

Restore bus number
Perform the second state estimator based on the result of the first estimator and PMU measurements ˆS x Calculate fitness of each individual by equations ( 29) and ( 30) The flow chart of the proposed optimal phasor measurement unit (PMU) placement (OPP) method in the compressed distribution network considering cost and state estimation accuracy.

Simulation Results and Discussion
The efficiency of the OPP method proposed in this paper is verified by the (Institute of Electrical and Electronics Engineers) IEEE standard bus systems.

Simulation Results and Discussion
The efficiency of the OPP method proposed in this paper is verified by the (Institute of Electrical and Electronics Engineers) IEEE standard bus systems.
The topology of IEEE 33 bus system, as shown in Figure 2, is used in this paper to carry out the analysis; 1-33 is the bus number.The SCADA measurements include the voltage amplitudes of buses 1, 2, 3, 6, 18, 23, 33; the current amplitudes of branches 3-4, 6-7, 2-19, 6-26; the injections at buses 5, 6, 13, 21; the power flows of branch 2-19, 23-24, 28-29.In the process of simulation, the power flow calculation results of the IEEE 33 bus system are selected as the truth values, and various measurements are generated by superimposing the corresponding random errors on the truth values.The standard deviations of measurement errors of voltage amplitude and current amplitude provided by the SCADA system are 1%.The standard deviation of the power measurement error is 2%.The standard deviations of measurement errors of voltage amplitude and branch current amplitude provided by PMU are 0.2%, the standard deviation of measurement errors of voltage phase angle provided by PMU is 0.5%, and the standard deviation of measurement errors of branch current phase angle provided by PMU is 0.5%.The standard deviation of the pseudo-measurement error is 20%.buses 5, 6, 13, 21; the power flows of branch 2-19, 23-24, 28-29.In the process of simulation, the power flow calculation results of the IEEE 33 bus system are selected as the truth values, and various measurements are generated by superimposing the corresponding random errors on the truth values.The standard deviations of measurement errors of voltage amplitude and current amplitude provided by the SCADA system are 1%.The standard deviation of the power measurement error is 2%.The standard deviations of measurement errors of voltage amplitude and branch current amplitude provided by PMU are 0.2%, the standard deviation of measurement errors of voltage phase angle provided by PMU is 0.5%, and the standard deviation of measurement errors of branch current phase angle provided by PMU is 0.5%.The standard deviation of the pseudo-measurement error is 20%.Before applying NSGA-II to solve the model, the network is compressed through the pre-configuration and existing measurements.Buses 5, 13, 21 are equivalent to "lesser ZIBs" with two neighbors, so buses 5, 13, 21 are merged with buses 4, 12, 22 into new buses 4', 12', 21' respectively.Bus 6 is equivalent to a "lesser ZIB" with three neighbors.Buses 1, 18, 21, 25, 33 are leaf buses, which are not placed with PMUs.The connected buses 2, 17, 20, 24, 32 are pre-configured with PMUs and merged with buses 1, 18, 21, 25, 33 into new buses 2', 17', 20', 24', 33' respectively.The compressed network and observable buses are shown in Figure 3. 23   The compressed system, as shown in Figure 3, is pre-configured with some PMUs.The buses pre-configured with PMUs are no longer involved in genetic operations, and their values of genes are 1.The specific parameters of applying NSGA-II to solve the model are as follows: The number of individuals in the parent population is 60.The number of individuals in the offspring population is 60.The number of iterations is 400.The number of buses in the system is 33.The number of branches is 32.The number of buses in the compressed system is 25.The number of variables in the placement vector is 25.The number of objectives is 2. The value of the initial front of each individual is 1.The number of intersections in the crossover operation is 2. The value of each subjective weight is 1.The number of rows in the chromosome matrix is 60, and the number of columns in the chromosome matrix is 29.The value of the penalty factor is 1000.Before applying NSGA-II to solve the model, the network is compressed through the pre-configuration and existing measurements.Buses 5, 13, 21 are equivalent to "lesser ZIBs" with two neighbors, so buses 5, 13, 21 are merged with buses 4, 12, 22 into new buses 4', 12', 21' respectively.Bus 6 is equivalent to a "lesser ZIB" with three neighbors.Buses 1,18,21,25,33   Before applying NSGA-II to solve the model, the network is compressed through the pre-configuration and existing measurements.Buses 5, 13, 21 are equivalent to "lesser ZIBs" with two neighbors, so buses 5, 13, 21 are merged with buses 4, 12, 22 into new buses 4', 12', 21' respectively.Bus 6 is equivalent to a "lesser ZIB" with three neighbors.Buses 1,18,21,25,33   The compressed system, as shown in Figure 3, is pre-configured with some PMUs.The buses pre-configured with PMUs are no longer involved in genetic operations, and their values of genes are 1.The specific parameters of applying NSGA-II to solve the model are as follows: The number of individuals in the parent population is 60.The number of individuals in the offspring population is 60.The number of iterations is 400.The number of buses in the system is 33.The number of branches is 32.The number of buses in the compressed system is 25.The number of variables in the placement vector is 25.The number of objectives is 2. The value of the initial front of each individual is 1.The number of intersections in the crossover operation is 2. The value of each subjective weight is 1.The number of rows in the chromosome matrix is 60, and the number of columns in the chromosome matrix is 29.The value of the penalty factor is 1000.The compressed system, as shown in Figure 3, is pre-configured with some PMUs.The buses pre-configured with PMUs are no longer involved in genetic operations, and their values of genes are 1.The specific parameters of applying NSGA-II to solve the model are as follows: The number of individuals in the parent population is 60.The number of individuals in the offspring population is 60.The number of iterations is 400.The number of buses in the system is 33.The number of branches is 32.The number of buses in the compressed system is 25.The number of variables in the placement vector is 25.The number of objectives is 2. The value of the initial front of each individual is 1.The number of intersections in the crossover operation is 2. The value of each subjective weight is 1.The number of rows in the chromosome matrix is 60, and the number of columns in the chromosome matrix is 29.The value of the penalty factor is 1000.
The Pareto set obtained by the optimization is shown in Figure 4.As shown in Figure 4, at least 10 PMUs are needed to ensure the observability of the network.The locations of PMUs are: buses 2, 8, 11, 15, 17, 20, 24, 26, 30, 32.Considering the cost and state estimation accuracy comprehensively, the compromise solution is selected by TOPSIS.The compromise solution needs 11 PMUs, and the locations are buses 2, 4, 8, 11, 15, 17, 20, 24, 26, 30, 32.Appl.Sci.2019, 9, x FOR PEER REVIEW 18 of 24 The Pareto set obtained by the optimization is shown in Figure 4.As shown in Figure 4, at least 10 PMUs are needed to ensure the observability of the network.The locations of PMUs are: buses 2, 8, 11, 15, 17, 20, 24, 26, 30, 32.Considering the cost and state estimation accuracy comprehensively, the compromise solution is selected by TOPSIS.The compromise solution needs 11 PMUs, and the locations are buses 2,4,8,11,15,17,20,24,26,30,32.In order to reflect the influence of network compression on OPP, this paper uses the single variable method to compare two cases, the result of the comparison is shown in Table 1.As shown in Table 1, on the one hand, network compression based on and existing measurements can effectively reduce the number of PMUs of the most economical solution and the compromised solution.On the other hand, under the same number of PMUs, the state estimation error of the solution is reduced by network compression.Note that: the most economical solution is the solution that has the minimum number of PMUs in the Pareto set; the compromised solution is the solution selected by the TOPSIS from the Pareto set.-4) In order to further analyze the influence of the compromised solution on the state estimation accuracy, 1000 state estimation simulations were carried out, the comparison of error in voltage estimation is shown Figure 5, the comparison of error in voltage phase angle is shown in Figure 6.In order to reflect the influence of network compression on OPP, this paper uses the single variable method to compare two cases, the result of the comparison is shown in Table 1.As shown in Table 1, on the one hand, network based on pre-configuration and existing measurements effectively the of PMUs of the most economical solution and the compromised solution.On the other hand, under the same number of PMUs, the state estimation error of the solution is reduced by network compression.Note that: the most economical solution is the solution that has the minimum number of PMUs in the Pareto set; the compromised solution is the solution selected by the TOPSIS from the Pareto set.In order to further analyze the influence of the compromised solution on the state estimation accuracy, 1000 state estimation simulations were carried out, the comparison of error in voltage amplitude estimation is shown in Figure 5, the comparison of error in voltage phase angle estimation is shown in Figure 6.
As shown in Figures 5 and 6, the compromised solution obtained by the OPP method proposed in this paper can not only reduce not only the voltage amplitude estimation error but also the voltage phase angle estimation error of each bus effectively.
In order to further analyze the effectiveness of searching the minimum number of PMUs required for the full observability, the simulation results of the proposed method in IEEE 33, 34 and 69 bus systems are compared with those of other methods in References [46][47][48], which are shown in Table 2.The topologies of the IEEE 34 bus system and IEEE 69 bus system are shown in Appendix A.  As shown in Figures 5 and 6, the compromised solution obtained by the OPP method proposed in this paper can not only reduce not only the voltage amplitude estimation error but also the voltage phase angle estimation error of each bus effectively.
In order to further analyze the effectiveness of searching the minimum number of PMUs required for the full observability, the simulation results of the proposed method in IEEE 33, 34 and 69 bus systems are compared with those of other methods in References [46][47][48], which are shown in Table 2.The topologies of the IEEE 34 bus system and IEEE 69 bus system are shown in Appendix A.    As shown in Figures 5 and 6, the compromised solution obtained by the OPP method proposed in this paper can not only reduce not only the voltage amplitude estimation error but also the voltage phase angle estimation error of each bus effectively.
In order to further analyze the effectiveness of searching the minimum number of PMUs required for the full observability, the simulation results of the proposed method in IEEE 33, 34 and 69 bus systems are compared with those of other methods in References [46][47][48], which are shown in Table 2.The topologies of the IEEE 34 bus system and IEEE 69 bus system are shown in Appendix A.  As shown in Table 2, compared with the methods in References [46][47][48], the OPP method proposed in this paper can provide the optimal solution with the minimum number of PMUs required for the full observability of IEEE 33, 34 and 69 bus systems respectively.Furthermore, the proposed OPP method can also obtain the compromise optimal scheme which comprehensively considers the cost and state estimation error.

Conclusions
Taking the cost and the state estimation accuracy into consideration, a multi objective method of optimal PMU placement is proposed in this paper.The following conclusions are obtained.The network compression based on the pre-configuration and existing measurements can not only reduce the number of PMUs of the most economical solution and the compromised solution, but also reduce the state estimation error under the same number of PMUs.Furthermore, the proposed OPP method can provide the most economical solution with the minimum number of PMUs required for the full observability and the compromised solution which considers the cost and state estimation accuracy comprehensively.Moreover, the compromised solution obtained by the proposed OPP method can effectively reduce the voltage amplitude estimation error and voltage phase angle estimation error.In addition, the most economical solution obtained by the proposed OPP method has the minimum number of PMUs required for the full observability than those provided by other methods.The expected future work will take the influence of network reconfiguration into consideration.

Conflicts of Interest:
The authors declare no conflict of interest.The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
. IEEE 34 bus system.
error covariances of the real and imaginary parts of branch current; i r wk and i x wk are the real and imaginary parts of branch current; I wk and δ wk are the amplitude and phase angle of branch current; σ 2 I wk and σ 2 δ wk −ϕi,se| n s.t.M ≥ b .

Merge 1 Whether
ZIBs and lesser ZIBs which have two neighbors Merge leaf bus with its neighbor Rearrange bus number Set Gen = 1 Output the solution Generate the initial population with Q individuals and choose it as the parent population Input parameters: A, SCADA measurements, ZIBs, iterations MaxGen, the number of individuals Q, and so on.Start Place PMU on bus whose neighbor is leaf bus Perform crossover operation and mutation operation to generate the offspring population Merge the offspring population with the parent population Choose the solution with the minimum number of PMUs as the optimal solution Calculate the rank and crowding distance of each individual Choose Q individuals according to the rank and crowding distance to generate the parent population Gen = Gen + Gen > MaxGen ？ Yes No Output the Pareto Set Whether the cost is limit ?Choose the solution from the Pareto set by equations (34)-(43) as the optimal

Figure 1 .
Figure 1.The flow chart of the proposed optimal phasor measurement unit (PMU) placement (OPP) method in the compressed distribution network considering cost and state estimation accuracy.

Figure 5 .
Figure 5.Comparison of error in voltage amplitude estimation.

Figure 6 .
Figure 6.Comparison of error in voltage phase angle estimation.

Figure 5 .
Figure 5.Comparison of error in voltage amplitude estimation.

Figure 5 .
Figure 5.Comparison of error in voltage amplitude estimation.

Figure 6 .
Figure 6.Comparison of error in voltage phase angle estimation.

Figure 6 .
Figure 6.Comparison of error in voltage phase angle estimation.
Author Contributions: X.K. and Y.W. conceived the idea of the study and conducted the research.Y.W. analyzed most of the data and wrote the initial draft of the paper.X.Y. and L.Y. contributed to finalizing this paper.Funding: This research was funded by the National Key Research and Development Program of China, grant number 2017YFB0902900 and 2017YFB0902902.
are leaf buses, which are not placed with PMUs.The connected buses 2, 17, 20, 24, 32 are pre-configured with PMUs and merged with buses 1, 18, 21, 25, 33 into new buses 2', 17', 20', 24', 33' respectively.The compressed network and observable buses are shown in Figure3., 13, 21; the power flows of branch 2-19, 23-24, 28-29.In the process of simulation, the power flow calculation results of the IEEE 33 bus system are selected as the truth values, and various measurements are generated by superimposing the corresponding random errors on the truth values.The standard deviations of measurement errors of voltage amplitude and current amplitude provided by the SCADA system are 1%.The standard deviation of the power measurement error is 2%.The standard deviations of measurement errors of voltage amplitude and branch current amplitude provided by PMU are 0.2%, the standard deviation of measurement errors of voltage phase angle provided by PMU is 0.5%, and the standard deviation of measurement errors of branch current phase angle provided by PMU is 0.5%.The standard deviation of the pseudo-measurement error is 20%.

Table 1 .
Comparison of two cases based on network compression.

Table 1 .
Comparison of two cases based on network compression.

Table 2 .
Comparison of optimal solutions in IEEE Standard Bus Systems.

Table 2 .
Comparison of optimal solutions in IEEE Standard Bus Systems.