Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Optimal Micro-PMU Placement for Improving State Estimation Accuracy via Mixed-integer Semidefinite Programming  PDF

  • Yang Peng
  • Zhi Wu (Member, IEEE)
  • Wei Gu (Senior Member, IEEE)
  • Suyang Zhou (Member, IEEE)
  • Pengxiang Liu (Graduate Student, Member)
School of Electrical Engineering, Southeast University, Nanjing, China

Updated:2023-03-25

DOI:10.35833/MPCE.2021.000615

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

Micro-phasor measurement units (μPMUs) with a micro-second resolution and milli-degree accuracy capability are expected to play an important role in improving the state estimation accuracy in the distribution network with increasing penetration of distributed generations. Therefore, this paper investigates the problem of how to place a limited number of μPMUs to improve the state estimation accuracy. Combined with pseudo-measurements and supervisory control and data acquisition (SCADA) measurements, an optimal μPMU placement model is proposed based on a two-step state estimation method. The E-optimal experimental criterion is utilized to measure the state estimation accuracy. The nonlinear optimization problem is transformed into a mixed-integer semidefinite programming (MISDP) problem, whose optimal solution can be obtained by using the improved Benders decomposition method. Simulations on several systems are carried out to evaluate the effective performance of the proposed model.

I. Introduction

PHASOR measurement units (PMUs) can provide real-time magnitude and phase angle information of voltage and current with high accuracy, from which many applications such as network observability, state estimation (SE), and safety protection and coordinated control can benefit [

1]. Currently, PMUs have been widely applied in transmission networks and the optimal PMU placement (OPP) problem for transmission networks has been well studied. However, the network structure of the distribution network is quite different from that of the transmission network. Most of distribution networks operate radially with a large number of nodes and low node-connectivity. Moreover, more and more distributed generations are connected to the distribution network, which results in rapid change of network state and much more complex faults [2]. However, on the one hand, some PMU placement methods in transmission networks are not well suited to distribution networks. On the other hand, the micro-phasor measurement unit (µPMU) with micro-second resolution and milli-degree accuracy capability has been developed, which is more suitable for complex distribution networks [3]. Thus, it is of great significance to study OPP considering the characteristics of distribution network and economic factor to facilitate the monitoring and controlling of the distribution network.

Most of existing OPP methods ensure complete network observability with a minimum number of installed PMUs or a minimum installation cost. The network observability is mainly divided into topological observability [

4] and numerical observability [5], where the former is usually achieved based on the decoupled measurement model and graph theory [6], and the latter is mostly achieved through the gain matrix of SE with full rank [7]. However, it needs to cover almost one third of the nodes to obtain complete network observability in distribution networks [8], which leads to large installation scale and high cost. Therefore, many OPP studies also take pseudo-measurements as well as traditional measurements from supervisory control and data acquisition (SCADA) system into consideration. This not only enhances the network observability, but also reduces the number of PMUs and the total cost for obtaining complete network observability. In recent literature, mathematical and heuristic optimization algorithms were implemented to solve OPP problem. Mathematical optimization algorithms include binary integer linear programming (BILP) [9], binary semidefinite programming (BSDP) [10], and recursive quadratic programming [11]. Interior-point method [12], branch and bound method [13], and Benders decomposition method [14] are commonly used to solve these programming problems. In [10], a BSDP approach was presented to solve the substation-oriented OPP problem considering the observability of transformer tap settings and the limited PMU channel capacity. Heuristic optimization algorithms include genetic algorithm (GA) [13], [15], simulated annealing algorithm (SAA) [16], binary particle swarm optimization (BPSO) [17], and evolutionary algorithm (EA) [18]. In [17], a BPSO-based methodology was proposed for the OPP to ensure topological observability while maximizing the measurement redundancy.

Furthermore, some other possible scenarios [

19] such as line outages [20], loss of PMUs [21], [22], controlled islanding [23], [24], and bounded observability propagation [25] were considered based on the network observability. A model for the optimal placement of contingency-constrained PMUs was presented in [20] by formulating conventional complete network observability first and then adding different contingency conditions in the network to the main model. In [23], an OPP model was proposed considering the controlled islanding in order that the network remains observable under both controlled islanding condition and normal operation condition. In [25], the optimal number and placement locations of PMUs were determined under new observability propagation rules to enhance measurement systems.

However, it is worth noting that studying the full rank of gain matrix or Jacobian matrix of SE can obtain network observability and ensure the executability of SE [

26], [27], but cannot improve the SE accuracy effectively. In recent years, some OPP methods were presented that take the improvement of SE accuracy as the optimization objective. Some researchers evaluated SE accuracy with scalar function of the SE covariance matrix [28] and proposed some OPP methods to improve the SE accuracy of distribution networks [29]-[37]. In [29], the mean squared error of SE subject to several observability constraints under a fixed number of PMUs was minimized. In [30], the concept of participation of various states in the error uncertainty was applied to identify the most important states that influence the estimated errors, and an OPP method was proposed without solving mathematical programing problems. In [31], the OPP problem was formulated as an optimal experiment design problem with E-, A- and D-optimal criteria to improve SE accuracy and the greedy approach was implemented to solve it, but the optimality of feasible solutions cannot be guaranteed. In [32], the D-optimal design was used to formulate OPP problem as a Boolean-convex optimization model based on fisher information matrix. In [33], the M-optimal experimental design technique was utilized to optimize the measurement allocation whereas the model cannot be solved efficiently. In [34], the OPP problems based on four optimal experimental criteria were formulated as semidefinite programming (SDP) problems by using convex relaxation, and a simple heuristic algorithm was proposed to obtain a feasible but not necessarily optimal binary solution. In [35], by choosing weighted least square (WLS) based SE under Gaussian hypothesis, a local optimization method following the convex relaxation was employed to solve the SDP problem. In [36], different distributions of noise measurement and robust estimators were considered, and the improved optimization solution method similar to [35] was utilized to obtain local optimal solution. In [37], an OPP method was proposed to minimize the mean squared error of SE measured by A-optimal criterion and feasible solutions were obtained by implementing Cholesky decomposition. Based on a two-step SE method, several bounds for the optimal solution to the PMU placement problem in distribution networks were compared in [38] by using greedy approach and convex relaxation, and the gap between a given suboptimal solution and the optimal solution can be checked with these bounds.

Therefore, on the basis of the existing studies, this paper proposes an optimal µPMU placement method to obtain the optimal solution, thus maximizing SE accuracy in distribution networks. A covariance matrix of SE error in rectangular coordinates is formulated based on mixed measurements and the two-step SE method [

39]. Then, an optimal µPMU placement model aimed at minimizing the worst error variance of SE is established based on E-optimal experimental criterion. After that, the optimization problem is transformed into a mixed-integer semidefinite programming (MISDP) problem which can be solved by using the improved Benders decomposition method.

The main contributions of this paper are listed as follows.

1) An optimal µPMU placement model is proposed to minimize the worst error variance of SE based on E-optimal experimental criterion for a given number of μPMUs, thus improving the SE accuracy. Pseudo-measurements, SCADA measurements, and zero injection buses (ZIBs) are collectively utilized for a two-step SE to reduce the number of required µPMUs.

2) The proposed nonlinear optimization model is transformed to an MISDP problem by introducing an auxiliary variable and the improved Benders decomposition method is used to obtain the optimization solution.

3) Global optimal solutions can be obtained by utilizing the improved Benders decomposition method compared with convex relaxation and greedy approach, and the solving time is much less than CUTSDP solver, conventional Benders decomposition method, and branch and bound method, especially when the number of μPMUs is large.

The rest of the paper is organized as follows. Section II introduces the proposed methodology for SE and an evaluation index of estimation accuracy. Section III establishes the optimal μPMU placement model and proposes the solution methodology. Section IV presents the results of case study. Finally, Section V draws the conclusions.

II. Proposed Methodology for SE and Evaluation Index of Estimation Accuracy

Among SE methods utilized in the OPP problem, the most common one is WLS method [

40]. A new two-step SE method considering mixed measurement and ZIB was proposed in [39] and was used to establish the OPP model in [38]. Compared with the traditional method, the two-step SE method avoids the recalculation of SE when the SCADA measurements change or multi-stage μPMU deployment is considered. This is because the SCADA measurements and prior installed μPMUs are processed in the second step after dealing with pseudo-measurements, whereas all the existing measurements are processed by the WLS-based method. Moreover, the proposed two-step SE method can effectively improve the computational efficiency of real-time SE while guaranteeing the accuracy of the results. Therefore, the two-step SE method is carried out to solve SE involved in optimal µPMU placement problem in this paper.

A. Two-step SE Method

For a three-phase network with N nodes (including Nsrc active nodes and Nload passive nodes), the brief calculation steps of the two-step SE method in rectangular coordinates are presented as follows. For more details, please refer to [

39].

In the first step, prior estimation is carried out by using pseudo-measurements, and the voltage formula at iteration k can be expressed as:

Vk=Bk-1R{Spsd}{Spsd}+V0 (1)

where R{*} and {*} are the real and imaginary parts of *, respectively; Spsd is the matrix of pseudo-measurements that are known beforehand; V0R6×Nload is the voltage matrix under zero loads in rectangular coordinates; and Bk-1 is the coefficient matrix at iteration k.

The iteration stops only when the gap between the calculated Spsd and the real Spsd is less than the given value. After the first step of SE, the estimation error covariance EpriorR(6Nload)×(6Nload) of VpriorC3×Nload can be obtained from (2).

Eprior=B0σpsddiag(R{Spsd})2({Spsd})2B0T (2)

where σpsd is the standard deviation of pseudo-measurements; and B0 is the initial coefficient matrix.

As analyzed in [

38] and [39], ZIBs lead to zero rows in Eprior thus making the matrix irreversible, and a linear transformation method is used to produce feasible solutions to this problem. But this method is unable to take µPMU placements at ZIBs into account because of dimensionality reduction of state variables. In this paper, zero value of the ZIB is replaced by a very small constant, i.e., 10-6, and a relatively small standard deviation, i.e., 0.01, is set for ZIBs. In this way, the µPMU placement of each node can be taken into consideration and the effect of approximate substitution method will be discussed in detail later.

In the second step, SCADA and µPMU measurements are utilized to carry out a linear post-estimate combined with Vprior and Eprior obtained in the first step. The post-estimate voltage Vpost can be expressed as:

Vpost=Vprior+KR{zP-CPVprior}{zP-CPVprior}zS-CSVprior (3)

where zPCNpm and zSRNsm are the vectors of Npm µPMU measurements and Nsm SCADA measurements, respectively; CP is the matrix mapping state voltages to µPMU measurements; CS is the matrix mapping state voltages to SCADA measurements; and K is the gain matrix obtained by minimizing the error covariance Epost, which can be expressed as:

K=EpriorCAT(CAEpriorCAT+EA)-1 (4)
CA=R{CP}-{CP}{CP}R{CP}VpriorCS(V)Vprior (5)
EA=diag(EP,ES) (6)

where EP and ES are the covariance matrices of the µPMU and SCADA measurement noises, respectively; and V is the state voltage vector.

The relationship between the voltage and µPMU measurement i at phase p of bus k is shown in (7), and the relationship between the voltage and SCADA measurement i at phase p of bus k is shown in (8).

(CPV)i=(CP)i,V=Vkp                               voltage measurement(Y)kp,V                       current measurement(Y)kp,mp(Vkp-Vmp)   current measurement of branch km (7)
(CS(V))i=|Vkp|                              voltage measurement|(Y)kp,V|                       current measurement|(Y)kp,mp(Vkp-Vmp)|   current measurement of branch km (8)

where Vkp is the voltage at phase p of bus k; (*)i denotes the ith element of vector *; (*)i, denotes the ith row of matrix*; and (*)kp,mp denotes the element of matrix * in the kpth row and mpth column; * denotes the magnitude of complex number *; and Y is the node admittance matrix of the system.

Finally, the first-order approximation of the estimation error covariance EpostR(6Nload)×(6Nload) can be obtained from (9).

EpostK(EA+CAEpriorCAT)KT+Eprior-KCAEprior-EpriorCATKT (9)

B. Evaluation Index of Estimation Accuracy

The error covariance matrix Epost is obtained through a two-step SE method in the previous subsection and indicates the estimation accuracy that the measurement system can achieve. The scalar-valued function f(Epost) is usually chosen as an indicator to evaluate the accuracy of SE in the estimation field [

37]. One of the functions is based on E-optimal design [28] and it can be expressed as:

fE=λmax(Epost) (10)

where λmax(*) is the maximum eigenvalue of * and its physical meaning is the worst error variance of an estimate. In this paper, the worst error variance, i.e., E-optimal design, is selected as the evaluation index of SE accuracy and the optimal µPMU placement model is established based on this strategy.

III. Optimal µPMU Placement Model and Solution Methodology

A. Optimal µPMU Placement Model

According to [

38], the deformation formula of Epost can be expressed as:

Epost=(Eprior-1+CA *EA-1CA)-1 (11)

where * denotes the complex conjugate transpose. Further, since the measurements at different nodes are independent of each other, CA *EA-1CA can be rewritten as:

CA *EA-1CA=isiim(CA')im,*(CA')im,(EA-1')im,im+iaiik(CA')ik,*(CA')ik,(EA-1')ik,ik (12)

where si and ai are the binary variables, si=1 means that a µPMU is installed at node i, and ai=1 means that other measurement device such as feeder terminal unit (FTU) is installed at node i, 0 otherwise; im is the serial number of the mth µPMU measurement of node i among all µPMU measurements; ik is the serial number of the kth SCADA measurement of node i among all SCADA measurements; and CA' and EA-1' are the values of CA and EA-1 when all measurements for all nodes are configured, respectively.

In this paper, we only focus on the optimal placement of µPMUs with unlimited number of channels, so the variable ai is known by default and only the variable si is optimized. Therefore, the latter formula in (12) related to ai can be replaced with a constant cSCADA and Epost can be expressed as (13) which is only related to the variable s=(si).

Epost(s)=Eprior-1+cSCADA+si(CA')im,*(CA')im,(EA-1')im,im-1 (13)

The objective of the optimal µPMU placement model in this paper is to improve the SE accuracy, and the constraint is the number of installed µPMUs. Therefore, combined with the E-optimal standard design in Section II-B as well as the SE error covariance matrix (13), the model can be written as:

minsfE=minsλmax(Epost(s))s.t. siNset    i       si{0,1}    i (14)

where Nset is the maximum number of µPMUs.

B. Solving Methodology

As shown in (13), Epost(s) in (14) involves the inversion of variable matrix, which makes (14) be a nonlinear optimization problem and hard to obtain the largest eigenvalue of Epost(s) directly. As we know, the maximum eigenvalue of a matrix equals to the inverse of the minimum eigenvalue of its inverse matrix. In this way, the original problem of minimizing the maximum eigenvalue of Epost(s) can be converted to a problem of how to maximize the minimum eigenvalue of Epost-1(s):

maxsfE=maxsλmin(Epost-1(s))s.t.  siNset    i        si{0,1}    i (15)

where λmin(*) is the minimum eigenvalue of *.

By introducing the auxiliary variable t [

34], the nonlinear E-optimal µPMU placement problem (15) can be rewritten as an MISDP problem (16), where variables si are all binary variables and variable t is a continuous variable.

min(-t)s.t. Epost-1(s)-tI2Nload0      siNset    i      si{0,1}    i (16)

where I2Nload is an identity matrix.

The proof of transformation is presented as follows.

Proof   Constraint Epost-1(s)-tI2Nload0 in (16) means that all eigenvalues ε of Epost-1(s)-tI2Nload are greater than or equal to 0, which means that when |Epost-1(s)-tI2Nload-εI2Nload|=0, ε0 always holds. In order to achieve this, the value of t needs to satisfy that the solution of the inequality t+ελmin(Epost-1(s)) with ε as a variable is greater than or equal to 0. In this case, when ε=εmin=0, t gets the maximum value which is equal to the minimum eigenvalue of Epost-1(s), i.e., λmin(Epost-1(s)). Therefore, (16) is equivalent to (15).

Considering the mathematical characteristics of the MISDP model (16), Benders decomposition is utilized in this paper to decompose the model into a simple mixed-integer linear programming (MILP) master problem and an SDP subproblem with continuous variables [

14]. The master problem determines the optimal location, and the subproblem optimizes the SE accuracy.

The subproblem is an SDP model where the discrete variables si are given, and the expression is as follows:

min(-t)s.t.  si=ξik:τik    i       Epost-1(s)-tI2Nload0 (17)

where ξik is the µPMU placement result of node i at iteration k; and τik is the dual variable corresponding to the constraint si=ξik at iteration k and it can be obtained when the subproblem is solved.

The master problem is an MILP model, and the expression is as follows:

min(α)s.t.  siNset    i        si{0,1}    i        α-tl+τil(si-ξil)    l=0,1,...,k (18)

where -tl is the objective function value of the subproblem with si=ξil at iteration l; and α is an auxiliary variable added in the solving process. During the process of solving the MISDP model based on the Benders decomposition for optimal µPMU placement, Benders cuts formed by -tl and τil play the role of connecting the master problem and the subproblem.

The master problem is formed by relaxing the constraints of the original problem, while the subproblem is a restricted version of the original problem. Therefore, the Benders decomposition method regards the objective function values of master problem and subproblem as lower and upper bounds of original problem, respectively. During iterations, the objective function value of the master problem increases with the Benders cuts provided by the subproblem. Finally, the method converges until the gap between objective function values of master problem and subproblem is less than a given value. However, it is worth noting that the new master problem must be solved from scratch when the Benders cut generated by the subproblem is added as a new constraint in each iteration. As the number of iterations increases, the more Benders cuts are added to (18), the longer it takes to solve the master problem.

Therefore, the improved Benders decomposition method based on lazy constraints is utilized to solve this problem in this paper [

41]. Lazy constraints are a group of inequalities specified by the user and are all necessary for the model. But these constraints are not added to the model initially. When a feasible solution is found and goes against the original problem, the corresponding lazy constraints are added to the model being solved currently. It can be observed that a model containing a large number of constraints can benefit from this method greatly.

The procedure of the improved Benders decomposition method using lazy constraints can be summarized as follows.

Step 1:   k=0. Start to solve the MILP master problem (18) without Benders cut. When a feasible solution ξk is found, go to Step 2.

Step 2:   solve the SDP subproblem (17) using feasible solution ξk to obtain the objective function value -tk of subproblem and the dual variable τk.

Step 3:   form the Benders cut according to the feasible solution ξk, the dual variable τk, and the objective function value -tk of subproblem. Add the Benders cut to the master problem as lazy constraint.

Step 4:   k=k+1. Continue to solve the master problem. If a feasible solution ξk is found, return to Step 2. If the convergence precision σ(10-4) of master problem is reached, output the optimal solution ξk.

With lazy constraints, there is no need to solve the master problem from scratch after each cut is added, which greatly accelerates the solving process. Therefore, using the improved Benders decomposition method to solve the MISDP problem ensures the solution efficiency.

In addition, it is known that the convergence of the improved Benders decomposition method for the MISDP problem is guaranteed as long as the envelope of function g(s1',s2',...,sNload') is convex, where g(s1',s2',...,sNload') is the function that relaxes all 0-1 variables into continuous variables and provides the optimal objective function value of problem (19) for given values of s1',s2',...,sNload' [

14].

g(s1',s2',...,sNload')=min(-t)s.t.  Epost-1(s')-tI2Nload0       si'Nset    i        si'[0,1]    i (19)

Since the SDP function g(s1',s2',...,sNload') is convex according to the definition of convex programming, the convergence of the proposed MISDP problem can be guaranteed, which ensures the optimal solution of discrete variables. Together with the fact that the SDP model (17) can ensure the global optimality of the solution of continuous variables, the proposed method can finally obtain σ-global optimal solution, i.e., the error between the obtained optimal solution and the true global optimal solution does not exceed σ.

To implement the proposed method, a complete flowchart is presented in Fig. 1.

Fig. 1  Flowchart of proposed method.

IV. Case Study

The proposed optimal µPMU placement model is evaluated in two balanced systems, namely the IEEE 33-bus and 123-bus systems [

42], which are shown in Figs. 2 and 3, respectively. And a practical distribution network with 446 nodes is also tested.

Fig. 2  Topology of IEEE 33-bus system.

Fig. 3  Topology of IEEE 123-bus system.

It is worth noting that balanced systems are considered for simplicity in this paper whereas the proposed method is also applicable to unbalanced distribution systems. The pseudo-measurements utilized in the test cases include active and reactive power injection of each node. Since these measurements are estimations rather than actual values, the corresponding noise is modeled with a relatively large standard deviation with σpsd=0.5 [

38], [39]. SCADA and µPMU measurements include voltage measurements, injected current measurements, and branch current measurements. The SCADA measurements can only provide magnitude information whereas the µPMU measurements can provide both magnitude and phase angle information. The standard deviation of the error is set to be 0.05 for SCADA measurements and 0.01 for µPMU measurements [43]. The simulations are performed in Python 3.7.9 [44] using an AMD R7-3700X CPU at 3.60 GHz with 32 GB of RAM. The master problem of the MISDP problem is solved by using Gurobi solver [45] and the Callback function in Gurobi is called to solve the subproblem and to add lazy constraints when a feasible solution of master problem is found. The SDP subproblem is solved by using Mosek solver [46] with CVXPY package [47].

A. µPMU Placement Results Under Different Combinations of ZIBs

In this subsection, approximate substitution method which is utilized to deal with ZIBs is compared with linear transformation method [

38] under different combinations of ZIBs in IEEE 33-bus system. Five ZIBs are randomly selected and added in order. The number of installed µPMUs Nset is set to be 4 and the SCADA measurements are not considered here. As shown in Table I, when 2, 4 and 5 ZIBs are considered, the approximate substitution method obtains larger objective values of the problem (15), which means the better optimization effect. For example, when nodes 6 and 11 are ZIBs, the approximate substitution method selects ZIB 11 as one of the µPMU installation nodes, which helps increase the objective value. In comparison, the linear transformation method neglects the µPMU placements at ZIBs because of dimensionality reduction. The same is true for the case with 4 ZIBs where ZIB 25 is selected to deploy the µPMU based on the approximate substitution method. The comparison indicates that it is better to take ZIBs into consideration during the µPMU placement.

TABLE I  μPMU Placement Results and Objective Values Using Approximate Substitution Method and Linear Transformation Method
No. of ZIBsApproximate substitution methodLinear transformation method [38]
µPMU placement resultObjective valueSolving time (s)µPMU placement resultObjective valueSolving time (s)
None 2, 6, 11, 25 62880 122 2, 6, 11, 25 62880 122
6 2, 7, 12, 26 75640 121 2, 7, 12, 26 75640 120
6, 11 2, 7, 11, 26 78000 128 2, 7, 12, 26 75780 112
6, 11, 14 1, 4, 9, 25 79710 116 1, 4, 9, 25 79710 108
6, 11, 14, 25 1, 5, 10, 25 90610 130 1, 4, 9, 26 82570 101
6, 11, 14, 25, 31 1, 5, 8, 13 114920 126 1, 3, 7, 9 97580 95

B. Sensitivity to SCADA Measurement Placement

The sensitivity of the µPMU placement results to the SCADA measurement placement is tested in IEEE 33-bus system. The number of nodes with SCADA measurements Ns ranges from 0 to 4. On the basis of the SCADA measurements in the previous round, one of the nodes configured with µPMUs from the previous round is selected to be equipped with SCADA measurements in the current round. ZIBs are not considered here. Figure 4 shows the µPMU placement results under different SCADA measurement placements using the proposed method when the number of installed µPMUs Nset is set to be 4.

Fig. 4  µPMU placement results when Nset=4 under different SCADA measurement placements using proposed method. (a) Ns=0. (b) Ns=1. (c) Ns=2. (d) Ns=3. (e) Ns=4.

It can be observed that the proposed method can take good account of the SCADA measurements in all cases, thus making the distribution of measurements in the system more uniform. Additionally, Fig. 5 shows the optimization results under different SCADA measurement placements and different numbers of µPMUs using the proposed method. A small maximum eigenvalue proves that the worst error variance is minor. The number of installed µPMUs Nset is limited to one quarter of the number of nodes considering economic factors, i.e., Nset ranges from 1 to 8. Figure 6 demonstrates that the more nodes are configured with SCADA measurements, the fewer µPMUs are required to achieve the same results. For instance, the result of Ns=4, Nset=3 is similar to those of Ns=3, Nset=4 and Ns=0, Nset=5. The result of Ns =2, Nset =3 is similar to that of Ns =1, Nset =4. However, as the number of µPMUs increases, the influence of SCADA measurements is gradually weakened. As shown in the Fig. 5, the gap between the five cases at Nset=8 is relatively smaller than that at Nset=4.

Fig. 5  Optimization results under different SCADA measurement placements and different numbers of µPMUs using proposed method.

Fig. 6  Optimization results using different methods in IEEE 33- and 123-bus systems. (a) IEEE 33-bus system. (b) IEEE 123-bus system.

C. Optimization Solution of MISDP Problem

Before tests, IEEE 33- and 123-bus systems are both configured with SCADA measurements and ZIBs randomly. Particularly, there are 11 SCADA measurements and 2 ZIBs in IEEE 33-bus system, and 43 SCADA measurements and 5 ZIBs in IEEE 123-bus system. Table II shows the types and quantities of SCADA measurements in two systems. The configuration of SCADA measurements and ZIBs above is kept unchanged for all the following tests.

TABLE II  Types and Quantities of SCADA Measurements in Two Systems
ItemIEEE 33-busIEEE 123-bus
No. of ZIBs 14, 30 14, 26, 36, 76, 110
No. of nodes with SCADA measurements 16, 19, 32 3, 12, 28, 42, 57, 69, 84, 89, 101, 122

Number of voltage

measurements

3 10
Number of injected current measurements 3 10
Number of branch current measurements 5 23

Firstly, the IEEE 33-bus system is tested in this case. Nset ranges from 1 to 8. Several methods are compared with the proposed method to solve the E-optimal µPMU placement problem.

1) Greedy Approach

The greedy approach in [

31] makes the optimal choice of µPMU placement at each round (totally Nset rounds). Specifically, in each round, the objective function value (14) is obtained for each candidate node configured with μPMU individually.

Then, the node with the minimum value is added to μPMU placement set and removed from candidate nodes. The approach terminates until the number of μPMUs reaches the upper limit.

2) Convex Relaxation Method

Convex relaxation method in [

34] is to relax the constraints si{0,1} in (16) to si[0,1], thus converting the MISDP problem to an SDP problem. After obtaining the optimal solution of SDP problem, the Nset largest variables of si are selected as μPMU placement nodes.

3) Enumeration Method

Enumeration method can list and compare all possible μPMU placement schemes, so it can ensure the optimality of final results. However, this method can only be utilized in small-scale systems considering long computation time.

The results of each method are shown in Fig. 6(a). The greedy approach performs the same as the proposed method for three out of eight times, but worse for the other times. All of the optimization results of the convex relaxation method are inferior to those of the other methods. Moreover, the results of the proposed method coincide with those of the enumeration method, which verifies its global optimality in IEEE 33-bus system.

Then, the IEEE 123-bus system is considered in this case. Nset is set as 1 to 30 and two methods, i.e., greedy approach and convex relaxation method, are compared with the proposed method. Figure 6(b) shows the results using different methods in IEEE 123-bus system. All of the optimization results of the convex relaxation method are inferior to the other methods obviously in 30 tests. Regarding greedy approach, it obtains the same optimal results as the proposed method in the first five rounds. But the gap gradually widens as the number of µPMUs increases. It can be observed that, similar to the IEEE 33-bus system, the proposed method outperforms the other methods in the IEEE 123-bus system, with better objective function values.

Additionally, Table III shows the optimal µPMU placement results in IEEE 33-bus system. It is noteworthy that the optimal sets of µPMU placement obtained by the proposed method do not satisfy the nested relations for different Nset. In other words, the optimal set of Nset µPMUs is not necessarily a subset of the optimal set of Nset+1 µPMUs. However, the number of the same µPMU placements between adjacent optimal sets increases significantly as Nset increases.

TABLE III  Optimal µPMU Placement Results in IEEE 33-bus System
Number of µPMUsOptimal µPMU placement result
1 5
2 2, 7
3 2, 5, 9
4 2, 6, 11, 25
5 1, 4, 7, 13, 28
6 1, 2, 5, 7, 13, 29
7 1, 2, 5, 7, 13, 25, 29
8 1, 2, 5, 7, 11, 14, 27, 29

Since the ultimate goal of optimal µPMU placement in this paper is to improve the SE accuracy, the effect of placement results on the SE accuracy is also investigated when Nset=5. Errors of voltage magnitude and phase angle in the IEEE 33-bus system are shown in Fig. 7. It can be observed that the proposed method yields a smaller relative error in both magnitude and phase angle of the estimated voltage than greedy approach and convex relaxation method, especially at nodes 4, 5, 16, 17, 18, 28, and 29.

Fig. 7  Errors of voltage magnitude and phase angle in IEEE 33-bus system. (a) Voltage magnitude. (b) Phase angle.

D. Applicability of Proposed Method to a Larger-scale System

A practical distribution network with 446 nodes is tested to discuss the applicability of the proposed method to practical power systems. A total number of 890 pseudo-measurements and 127 SCADA measurements are considered. Nset is set to be 1, 4, 7, , 100, respectively. Greedy approach and convex relaxation method are compared with the proposed method. It can be observed from Fig. 8 that the results of the convex relaxation method are inferior to the other methods and the results of the greedy approach are between those of the other two methods. The better optimization results of the proposed method demonstrate that the proposed method is available in a relatively large-scale power system.

Fig. 8  Optimization results using different methods in a practical distribution network.

E. Benefits of Improved Benders Decomposition Method

Finally, the E-optimal µPMU placement problem modeled as an MISDP problem is also solved by utilizing the following three methods: ① the conventional Benders decomposition (Method 1) introduced in Section III; ② the CUTSDP solver [

48] (Method 2); and ③ branch and bound method [49] (Method 3).

Method 2 is provided by a MATLAB toolbox named YALMIP [

50] for mixed-integer semidefinite programs, implementing an outer approximation approach. It is based on an iterative process where it relaxes the conic constraints to linear constraints, solves an MILP problem, and adds violated linear cuts to the model. Method 3 splits the initial MISDP problem into smaller subproblems, which are then solved in some specified order. Method 3 is implemented based on the BNB solver [51] for MISDP.

Table IV compares the solving time of the MISDP problem using Methods 1, 2, 3, and the proposed method in IEEE 33-bus system. It can be observed that the objective function values of the proposed method are consistent with the other methods. When Nset=1, 2, 3, Method 2 takes a similar amount of time as the proposed method whereas Method 1 takes nearly four times longer and Method 3 takes the longest time among all methods. When Nset exceeds 3, the solving time for Methods 1, 2, and 3 increases sharply. Particularly, Methods 1, 2, and 3 take 37.8 hours, 79.4 hours, and 5.2 hours to optimize the placement of 6 µPMUs, respectively, whereas the proposed method takes only 12 min. Moreover, the convergence diagrams of the four methods are presented in Fig. 9 when the number of µPMUs is 4. It can be observed that the convergence time of the proposed method is much less than the others. This is mainly because lazy constraints of the proposed method avoid the solving of master problem from scratch after each cut is added, while Method 1 and Method 2 solve the respective master problem repeatedly. In addition, since Method 3 is a deterministic method of search and iteration which belongs to partial enumeration, the solving time is longer when it is used individually to solve a real power system with many nodes.

TABLE IV  Comparison of Solution Information in IEEE 33-bus System
NsetProposed methodMethod 1Method 2Method 3
Solving time (s)Objective valueSolving time (s)Objective valueSolving time (s)Objective valueSolving time (s)Objective value
1 2 3800 8 3800 1 3800 47 3800
2 18 9000 86 9000 11 9000 255 9000
3 33 54300 183 54300 51 54300 637 54300
4 130 65800 1104 65800 1654 65800 2883 65800
5 399 123100 53180 123100 17860 123100 7290 123100
6 724 192500 285700 192500 136300 192500 18500 192500
7 809 211800 48510 211800
8 2775 323200 70200 323200

Fig. 9  Convergence diagram of four methods in IEEE 33-bus system.

Moreover, the solving time of the proposed method in the case of 30 µPMUs is 3.1 hours and 11.5 hours for the IEEE 123-bus system and the practical distribution network, respectively. While Methods 1, 2, and 3 cannot complete the solution of the same OPP problem within 24 hours. In summary, the proposed method shows higher solution efficiency than Methods 1, 2, and 3. It proves that the proposed method can deal with discrete variables and continuous variables efficiently and can obtain the global optimal solution at the same time.

V. Conclusion

This paper presents an optimal µPMU placement method in distribution networks. A new model for optimal µPMU placement is established based on a two-step SE method, which can not only improve the SE accuracy, but also take into account the impact of SCADA measurements and ZIBs. The nonlinear optimization problem is transformed into an MISDP problem and an improved Benders decomposition method is utilized to solve it to obtain global optimal solution. Simulations on several test cases verify the feasibility, optimality, and high efficiency of the proposed method by contrast with other methods.

Future work could include studying the optimal µPMU placement problem by combining SE with other applications such as fault location to meet multi-application requirements of distribution networks.

References

1

A. G. Phadke and T. Bi, “Phasor measurement units, WAMS, and their applications in protection and control of power systems,” Journal of Modern Power Systems and Clean Energy, vol. 6, no. 4, pp. 619-629, Jul. 2018. [Baidu Scholar] 

2

J. Xu, B. Xie, S. Liao et al., “Load shedding and restoration for intentional island with renewable distributed generation,” Journal of Modern Power Systems and Clean Energy, vol. 9, no. 3, pp. 612-624, May 2021. [Baidu Scholar] 

3

E. Dusabimana and S. G. Yoon. “A survey on the micro-phasor measurement unit in distribution networks,” Electronics, vol. 9, no. 2, p. 305, Feb. 2020. [Baidu Scholar] 

4

K. Chauhan and R. Sodhi, “Placement of distribution-level phasor measurements for topological observability and monitoring of active distribution networks,” IEEE Transactions on Instrumentation and Measurement, vol. 69, no. 6, pp. 3451-3460, Jun. 2020. [Baidu Scholar] 

5

F. Wu and A. Monticelli, “Network observability: theory,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-104, no. 5, pp. 1042-1048, May 1985. [Baidu Scholar] 

6

A. Ahmadi, Y. Alinejad-Beromi, and M. Moradi, “Optimal PMU placement for power system observability using binary particle swarm optimization and considering measurement redundancy,” Expert Systems with Applications, vol. 38, no. 6. pp. 7263-7269, Jun. 2011. [Baidu Scholar] 

7

N. M. Manousakis and G. N. Korres, “Optimal PMU placement for numerical observability considering fixed channel capacity – a semidefinite programming approach,” IEEE Transactions on Power Systems, vol. 31, no. 4, pp. 3328-3329, Jul. 2016. [Baidu Scholar] 

8

Z. Wu, X. Du, W. Gu et al., “Optimal PMU placement considering load loss and relaying in distribution networks,” IEEE Access, vol. 6, pp. 33645-33653, May 2018. [Baidu Scholar] 

9

V. S. Devendran, J. Jasni, M. A. Mohd Radzi et al., “Optimal placement of PMU for complete observability of power system considering zero injection and islanding condition,” in Proceedings of 2020 IEEE International Conference on Power and Energy (PECon), Penang, Malaysia, Dec. 2020, pp. 107-112. [Baidu Scholar] 

10

N. M. Manousakis and G. N. Korres, “Optimal PMU arrangement considering limited channel capacity and transformer tap settings,” IET Generation, Transmission & Distribution, vol. 14, no. 24, pp. 5984-5991, Nov. 2020. [Baidu Scholar] 

11

N. P. Theodorakatos, M. Lytras, and R. Babu, “Towards smart energy grids: a box-constrained nonlinear underdetermined model for power system observability using recursive quadratic programming,” Energies, vol. 13, p. 1724, Apr. 2020. [Baidu Scholar] 

12

N. P. Theodorakatos, “A nonlinear well-determined model for power system observability using interior-point methods,” Measurement, vol. 152, p. 107305, Nov. 2019. [Baidu Scholar] 

13

N. P. Theodorakatos, “Optimal phasor measurement unit placement for numerical observability using branch-and-bound and a binary-coded genetic algorithm,” Electric Power Components and Systems, vol. 47, no. 4-5, pp. 357-371, May 2019. [Baidu Scholar] 

14

A. J. Conejo, E. Castillo, R. Mínguez et al., Decomposition Techniques in Mathematical Programming. Berlin: Springer, 2006, pp. 243-270. [Baidu Scholar] 

15

H. H. Muller and C. A. Castro, “Genetic algorithm-based phasor measurement unit placement method considering observability and security criteria,” IET Generation, Transmission & Distribution, vol. 10, no. 1, pp. 270-280, Jan. 2016. [Baidu Scholar] 

16

P. Yuan, Q. Ai, and Y. Zhao, “Research on multi-objective optimal PMU placement based on error analysis theory and improved GASA,” Proceedings of the CSEE, vol. 34, no. 13, pp. 2178-2187, May 2014. [Baidu Scholar] 

17

N. H. A. Rahman and A. F. Zobaa, “Integrated mutation strategy with modified binary PSO algorithm for optimal PMUs placement,” IEEE Transactions on Industrial Informatics, vol. 13, no. 6, pp. 3124-3133, Dec. 2017. [Baidu Scholar] 

18

S. Prasad and D. M. V. Kumar, “Trade-offs in PMU and IED deployment for active distribution state estimation using multi-objective evolutionary algorithm,” IEEE Transactions on Instrumentation and Measurement, vol. 67, no. 6, pp. 1298-1307, Jun. 2018. [Baidu Scholar] 

19

M. U. Usman and M. O. Faruque, “Applications of synchrophasor technologies in power systems,” Journal of Modern Power Systems and Clean Energy, vol. 7, no. 2, pp. 211-226, Mar. 2019. [Baidu Scholar] 

20

F. Aminifar, A. Khodaei, M. Fotuhi-Firuzabad et al., “Contingency-constrained PMU placement in power networks,” IEEE Transactions on Power Systems, vol. 25, no. 1, pp. 516-523, Feb. 2010. [Baidu Scholar] 

21

R. Emami and A. Abur, “Robust measurement design by placing synchronized phasor measurements on network branches,” IEEE Transactions on Power Systems, vol. 25, no. 1, pp. 38-43, Feb. 2010. [Baidu Scholar] 

22

X. Kong, X. Yuan, Y. Wang et al., “Research on optimal D-PMU placement technology to improve the observability of smart distribution networks,” Energies, vol. 12, no. 22, pp. 1-23, Nov. 2019. [Baidu Scholar] 

23

L. Huang, Y. Sun, J. Xu et al., “Optimal PMU placement considering controlled islanding of power system,” IEEE Transactions on Power Systems, vol. 29, no. 2, pp. 742-755, Mar. 2014. [Baidu Scholar] 

24

S. Akhlaghi, N. Zhou, and N. Wu, “PMU placement for state estimation considering measurement redundancy and controlled islanding,” in Proceedings of 2016 IEEE PES General Meeting, Boston, USA, Jul. 2016, pp. 1-5 [Baidu Scholar] 

25

X. Guo, C. Liao, and C. Chu, “Enhanced optimal PMU placements with limited observability propagations,” IEEE Access, vol. 8, pp. 22515- 22524, Jan. 2020. [Baidu Scholar] 

26

R. Sodhi, S. C. Srivastava, and S. N. Singh, “Optimal PMU placement method for complete topological and numerical observability of power system,” Electric Power Systems Research, vol. 80, no. 9, pp. 1154-1159, Sept. 2010. [Baidu Scholar] 

27

G. N. Korres, N. M. Manousakis, T. C. Xygkis et al., “Optimal phasor measurement unit placement for numerical observability in the presence of conventional measurements using semi-definite programming,” IET Generation, Transmission & Distribution, vol. 9, no. 15, pp. 2427-2436, Nov. 2015. [Baidu Scholar] 

28

F. Pukelsheim, Optimal Design of Experiments. New York: Wiley, 1993, pp. 210-246. [Baidu Scholar] 

29

Y. Shi, H. D. Tuan, T. Q. Duong et al., “PMU placement optimization for efficient state estimation in smart grid,” IEEE Journal on Selected Areas in Communications, vol. 38, no. 1, pp. 71-83, Jan. 2020. [Baidu Scholar] 

30

E. Caro, R. Singh, B.C. Pal et al., “Participation factor approach for phasor measurement unit placement in power system state estimation,” IET Generation, Transmission & Distribution, vol. 6, no. 9, pp. 922-929, Sept. 2012. [Baidu Scholar] 

31

Q. Li, R. Negi, and M. D. Ilic, “Phasor measurement units placement for power system state estimation: a greedy approach,” in Proceedings of 2011 IEEE PES General Meeting, Detroit, USA, Jul. 2011, pp. 1-8. [Baidu Scholar] 

32

T. C. Xygkis, G. N. Korres, and N. M. Manousakis, “Fisher information-based meter placement in distribution grids via the D-optimal experimental design,” IEEE Transactions on Smart Grid, vol. 9, no. 2, pp. 1452-1461, Mar. 2018. [Baidu Scholar] 

33

T. C. Xygkis and G. N. Korres, “Optimized measurement allocation for power distribution systems using mixed integer SDP,” IEEE Transactions on Instrumentation and Measurement, vol. 66, no. 11, pp. 2967-2976, Nov. 2017 [Baidu Scholar] 

34

V. Kekatos, G. B. Giannakis, and B. Wollenberg, “Optimal placement of phasor measurement units via convex relaxation,” IEEE Transactions on Power Systems, vol. 27, no. 3, pp. 1521-1530, Aug. 2012 [Baidu Scholar] 

35

Y. Yao, X. Liu, and Z. Li, “Robust measurement placement for distribution system state estimation,” IEEE Transactions on Sustainable Energy, vol. 10, no. 1, pp. 364-374, Jan. 2019. [Baidu Scholar] 

36

T. Chen, Y. Cao, X. Chen et al., “Optimal PMU placement approach for power systems considering non-Gaussian measurement noise statistics,” International Journal of Electrical Power & Energy Systems, vol. 126, p. 106577, Mar. 2021. [Baidu Scholar] 

37

J. Tian, D. Liang, L. Ge et al., “Placement of micro-phasor measurement units in distribution systems for highly accurate state perception,” Power System Technology, vol. 43, no. 7, pp. 2235-2242, Jul. 2019. [Baidu Scholar] 

38

M. Picallo, A. Anta, and B. de Schutter, “Comparison of bounds for optimal PMU placement for state estimation in distribution grids,” IEEE Transactions on Power Systems, vol. 34, no. 6, pp. 4837-4846, Nov. 2019. [Baidu Scholar] 

39

M. Picallo, A. Anta, A. Panosyan et al., “A two-step distribution system state estimator with grid constraints and mixed measurements,” in Proceedings of 2018 Power Systems Computation Conference (PSCC), Dublin, Ireland, Jun. 2018, pp. 1-7. [Baidu Scholar] 

40

A. Abur and A. Gómez Expósito, Power System State Estimation: Theory and Implementation. Boca Raton: CRC, 2004, pp. 24-51. [Baidu Scholar] 

41

S. Lin, G. J. Lim, and J. F. Bard, “Benders decomposition and an IP-based heuristic for selecting IMRT treatment beam angles,” European Journal of Operational Research, vol. 251, no. 3, pp. 715-726, Jan. 2016. [Baidu Scholar] 

42

IEEE PES AMPS DSAS Test Feeder Working Group. (2021, May). IEEE 33-bus test system and IEEE 123-bus test system. [Online]. Available: http://sites.ieee.org/pes-testfeeders/resources/ [Baidu Scholar] 

43

C. Zhang, Y. Jia, Z. Xu et al., “Optimal PMU placement considering state estimation uncertainty and voltage controllability,” IET Generation, Transmission & Distribution, vol. 11, no. 18, pp. 4465-4475, Sept. 2017. [Baidu Scholar] 

44

Python. (2021, Jun.). The Python tutorial. [Online]. Available: https://docs.python.org/release/3.6.9/tutorial/index.html [Baidu Scholar] 

45

Gurobi. (2021, Jun.). Gurobi optimizer reference manual. [Online]. Available: https://www.gurobi.com/documentation/9.1/refman/index.html [Baidu Scholar] 

46

MOSEK. (2021, Jul.). MOSEK optimizer API for Python. [Online]. Available: http://www.mosek.com [Baidu Scholar] 

47

S. Diamond and S. Boyd, “CVXPY: a Python-embedded modeling language for convex optimization,” Journal of Machine Learning Research, vol. 17, no. 1, pp. 2909-2913, Apr. 2016. [Baidu Scholar] 

48

YALMIP. (2021, Jun.). The CUTSDP solver explained. [Online]. Available: https://yalmip.github.io/The-cutsdp-solver [Baidu Scholar] 

49

C. Rowe and J. Maciejowski, “An efficient algorithm for mixed integer semidefinite optimization,” in Proceedings of IEEE American Control Conference, Denver, USA, Jun. 2003, pp. 4730-4735. [Baidu Scholar] 

50

J. Lofberg, “YALMIP: a toolbox for modeling and optimization in MATLAB,” in Proceedings of 2004 IEEE International Conference on Robotics and Automation, Taipei, China, Sept. 2004, pp. 284-289. [Baidu Scholar] 

51

YALMIP. (2021, May). The BNB solver explained. [Online]. Available: https://yalmip.github.io/solver/bnb/ [Baidu Scholar]