An Optimization-Based Robust Dynamic State Estimation for Power Systems with Synchronized Phasor Measurement Units, Involving Disturbance Rejection

,


Introduction
Te operating situation of a power system at any time is considered the state of that power system. In other words, the state of a power system refects the behavior of that system. Te power system can be in any of the three states, normal state, emergency state, and restorative state. Te ultimate goal of the system operator is to keep the power system in normal condition. Achieving this goal requires system security analysis [1].
State estimation is a mathematical process that estimates the state variables (magnitude and angle of voltage of busbars) through system data and meter information. Utilizing state estimation in an energy management system (EMS) is important to accomplish various planning and control tasks [1]. Tere are mainly two types of state estimators; static (SSE) and dynamic (DSE). SSE is not able to be implemented in short time intervals, while DSE (or forecasting-aided state estimation) can estimate the state of the power system in short time steps, and in each time step, it uses two sets of measured and predicted data [2]. DSE is considered more practical for performing EMS tasks. Te DSE can predict the state of the power system gradually over short periods. In dynamic state estimation, unlike static state estimation, in which only a set of measurements is used to estimate the network state variables, two sets of measurements and predictions are used to achieve the goal of the state estimation process. Measurements are used to calculate predictions, and there are several methods for this purpose. One of the most common methods is the Holt method or the double exponential smoothing method.
State estimation utilizes available measurements. CVTs, CTs, relays, and PMUs are important meters used in power systems. PTs and CTs measure high-level voltages and currents and convert them into operational safe surface signals. Relays and PMUs are powered by CTs and PTs. Te measurements must be synchronized chronologically to obtain an accurate picture of the power system [3].
Te supervisory control and data acquisition (SCADA) system is a type of control system that collects data from various locations remotely, and through the collected data, monitors the power system. Data in digital or analog format are obtained from various meters located at substations. Te data are transmitted by meters to remote terminal units (RTUs) in the substation. RTU transmits data (including voltage, current, and circuit breaker status.) from the substation to computers in the control room via a telecommunications network. Te information received from the RTU is processed by a state estimator [4]. Te host computer in SCADA, which is located in the control room, receives all measurements through one of the available types of communication equipment (such as fber optics, satellite, and microwaves) and transmits them to the EMS via the local area network (LAN) [1].
Several capabilities are assumed for the new generation of power system monitoring schemes using the wide area monitoring system (WAMS) with the intensifcation of the infuence of phasor measurement units (PMUs). Tese capabilities include more accurate measurements from PMUs, higher sampling rates, temporal synchronization of data, and linearity between measurement values and state variables. Tus, the efect of considering PMUs in the state estimation process should be analyzed. In recent years, several methods have been proposed to improve and develop the performance of this issue [5]. WAMS complements the various tasks performed in the feld of data acquisition by sections such as protection relays, fault recorders, and SCADA, and can install PMUs and take advantage of the global positioning system (GPS).
Recent research has taken diferent approaches to the subject of state estimation of PMUs, which fall into several main categories of modifcation of state estimation methods (including robust state estimation, timing correlation, utilizing neural networks, state estimation optimization, and visibility analysis), according to the type of system (including decentralized wide area system and multiarea system), and with taking into account other elements of the power system (including HVDC communication lines and FACTS devices).
References [6][7][8][9][10][11][12][13][14][15] focus on the robustness of state estimation. It is worth noting that while most related studies are inherently based on the assumption that state estimation noise is known and has a Gaussian behavior [16], some others, such as [17][18][19][20], consider the non-Gaussian noise and ignore the possibility of occurring disturbances. In [7], system states are estimated using a robust dynamic state estimator based on PMUs. Tis estimator combines the existing SCADA data with accurate data obtained from PMUs using a method called fusion measurements based time-variant state transition matrix updating method (FSTMU), to instantly monitor the power system under various conditions; thus, the accuracy of the DSE process is increased. Te idea of [8] is to provide an adaptive weighting function to adjust the weight of each of the measurements generated by the PMU based on the spatial distance of the site of large unwanted disturbances to that measurement. Similar to [8], the proposed methods in [2,21,22] have the ability of disturbance rejection; however, their proposed methods ignore the concept of robustness against outliers or ignore the variety of diferent kinds of network sudden changes (disturbances). Te method used in [6] is a generalization of a set of several robust state estimation methods including quadratic constant (QC), quadratic linear (QL), square-root (SR), and multiple-segment (MS) methods. By combining these methods with the moving horizon estimation (MHE) method, a method called reweighted moving horizon estimation (RMHE) is introduced. Reference [9] performs the state estimation process in terms of network imbalance using modal decomposition. Reference [10] uses an imperial competitive algorithm (ICA) to make the dynamic state estimation process robust. Compared to other methods based on statistics or the Kalman flter, this method has some advantages, such as the lack of dependence on system linearity and Gaussian noise, resistance to falsifed data due to the assumptions made in the formation of imperialists in the ICA method, and due to the use of predictability in the proposed method. In [11], similar to [6,13,14], the noise of measurements is assumed to be non-Gaussian. In this reference, a two-stage state estimation process is presented, the frst stage of which uses a generalized maximum probability estimator (GM estimator). Tis GM estimator uses traditional SCADA measurements. In the second stage, the results of the frst stage are combined with PMU measurements to achieve a linear and robust estimate.
As other available methods of modifying the state estimation process, [2,[23][24][25][26][27][28][29][30][31] utilize neural networks and state estimation optimization. For example, [31] utilizes the stability and robustness of football game optimization (FGO) as a swarm-intelligence-based optimization technique in state estimation problems. In the application of estimating the dynamic state of power systems, since the phasor measurement units do not extend to all branches of the power system, the use of forecast data as the second set of DSE process inputs along with measurements increases the redundancy of input data to the process. Reference [2] provides a mixed-integer program (MIP) formulation for dynamic state estimation that can simultaneously discard predicted invalid values when sudden changes in the system state are detected. As another method of modifying the state estimation process, [32,33] analyze visibility. Reference [34] estimates the power system state in the presence of traditional HVDC communication lines with operation under diferent control modes, and [35] estimates the state in the presence of FACTS and devices of unifed power fow control (UPFC) type. [36][37][38][39][40] focus on state estimation according to the type of system (including decentralized wide-area systems and multiarea systems). For example, [38] provides a multiarea state estimation (MASE) approach. Tis approach is suitable for decentralized interconnected networks and uses the common external equivalent technique in terms of PMU measurements. A comprehensive review of DSE motivations, defnitions, and methods is available in [41]. It is worth mentioning that a wide range of research (e.g., see [42]) has also focused on enhancing the state estimation process regarding cyber security issues. Moreover, some references such as [43][44][45][46] focus on important related topics such as fault location in the distribution network.
As a summary of the topic, the image obtained from the power system must be improved to control the network quantities in the event of unwanted changes or disturbances in the grid. To achieve this, due to the high rate of change of data and the high volume of exchange information, the use of phasor measurement units can play an efective role. On the other hand, due to the inevitability of equipment measurement errors, to provide a favorable estimate of system states using the measured data, it is necessary to use a state estimator that can eliminate or reduce the adverse efect of such defective or bad data on the fnal output of the process which is called the robust state estimator. Te lack of enough related research at the transmission level of power systems considering the simultaneous signifcant efects of outliers, non-Gaussian noise of measurement, and disturbances, led the authors to thoroughly analyze this important issue. In this paper, in addition to using a dynamic state estimator with the ability to reduce the efect of disturbances on the network, the robustness of the state estimation process is also considered against the defective data afecting the existing measurements. Moreover, unlike the conventional repetitive algorithms and methods, an optimizationbased and two-stage method is proposed, which simultaneously increases the redundancy of inputs of the process and utilizes the positive features of the method of dynamic state estimation and robust state estimation. Tis proposed method can merge the traditional measurements and the ones obtained from PMUs. Moreover, it hastens the execution of the DSE and can simultaneously override the estimation process predictions when the estimated states exceed the estimation confdence bound. By implementing an explanatory three-bus network and the IEEE standard 14bus and 57-bus networks using the MATPOWER tool, simulation results are thoroughly analyzed.
Te rest of this paper is presented as follows. Section 2 provides an overview of the DSE and RSE. Te proposed RDSE model and formulation are presented in Section 3. Section 4 discusses numerical examples, and the fnal points are stated in Section 5.

Robust Dynamic State Estimation
2.1. Dynamic State Estimation. Dynamic state estimation (DSE) determines the current state of the system based on the values of the previous step state. Tus, unlike static state estimation, which belongs to the maximum probability estimation, dynamic state estimation is an example of Bayesian estimation in which recursive algorithms such as the Kalman flter are used to estimate the system states. Te accuracy of the DSE depends on the sampling rate of the measurements. Te low sampling rate of the traditional SCADA system does not allow the dynamic state of the system to be recorded, and therefore, the DSE cannot be implemented based on those measurements. Te introduction of the PMU with a high sampling rate between 30-60 Hz made it possible to implement the Kalman flter algorithm to implement DSE [4].
Te power system follows a quasistatic behavior under normal operating conditions; the daily profle is subject to small random fuctuations. Under these conditions, the power system state changes slowly, and the DSE performance is acceptable. However, during various disturbances (such as load fuctuations, generator outages, and network switchings), the state variables corresponding to the busbars that are close to the area afected by the disturbance undergo large transient changes. Under these conditions, the difference between the measurement values and the data derived from them, such as the predictions, will be large, and consequently, the adverse efect of the invalid predictions on the output values of the estimation state will reduce the performance of the dynamic estimation process.
Te method mentioned in this research uses measurements to calculate the forecast data set and then enters the two measurement and prediction sets as input to the state estimation problem. Tis process is based on an optimization method. In this optimization method, several measurement values are related to a state variable corresponding to each of the system busbars, based on the adjacent buses connected to it. Moreover, the degree of validity for each of the predictions and measurements in the objective function of the optimization problem is determined based on their standard deviations. Furthermore, the potential uncertainties in the problem are taken into account. As a result, it is possible to reduce the adverse efect of potential disturbances to the network on the accuracy of predictions, which consequently, improves the performance of the estimation process.
As mentioned in the DSE method, the problem inputs consist of two sets of data; measurement and prediction. In the method mentioned in this section, the measurements themselves are divided into two categories; direct measurements and indirect measurements. Interpretation of direct measurements refers to the phase and voltage values measured in PMU-equipped buses. Indirect measurements, on the other hand, refer to the measured values of the voltage phase (magnitude and voltage angle) in PMU-free buses.
Indirect measurements are calculated using a specifc formulation called GUM as an abbreviation. In this formulation, in addition to calculating the size and voltage phase of the bus without PMU (indirect measurements), the amount of standard deviation of each is also determined. In determining the standard deviation, there are uncertainties in the problem, including uncertainties corresponding to the measurements obtained from PMUs (including the magnitude and phase of voltage and magnitude and phase of current) and uncertainties corresponding to the transmission line parameters (including the magnitude and phase of the series impedance of the lines and the magnitude and International Transactions on Electrical Energy Systems phase of parallel admittance of the transmission lines) which are taken into account.
In the DSE method, the most common method of calculating predictions is used, which is the Holt linear exponential smoothing method [47]. Te Holt method used in this section is a conventional representation of the Holt method used in several kinds of research.

Robust State Estimation.
One of the purposes of using the state estimator is to detect and eliminate potential errors in measurements, network models, or parameters. If the estimated state of the power system is not afected by deviations from the true values of several measurements, then the corresponding estimator is considered robust. Te robustness of the state estimation process can often be implemented as the computations become more complex. Robust estimators must provide good output against different types of outliers [1].
Noise (error) related to the state estimation of the power system and examined in most related research works, is assumed to be Gaussian (normal); however, when the network is exposed to transient data in steady-state measurements, measurement equipment error, human error, or nonlinear factors, the measurement noise follows a non-Gaussian function. In such cases, defective data or outliers increase the likelihood of the improper performance of the estimation process. To overcome the inappropriate efect of the mentioned adverse data on the output of the state estimation process, various methods such as robust state estimation (RSE) have been proposed [6]. Te robustness of the estimation process means the ability to reduce or eliminate the impact of the process outliers. It is necessary to emphasize that bad data detection is another concept that is not in the purpose of this paper.
Considering recent research, it is necessary to propose a robust dynamic state estimation (RDSE) process, which simultaneously reduces or eliminates the efect of outliers due to measurement errors, while invalid predictions due to disturbances in the network are withdrawn.
It is worth mentioning that measurement noise (error) and outliers are two diferent concepts; i.e., the frst case shows the diference between the measured values for a state variable and its true value, while the second case represents the undesirable data that are considered as measurements. For example, when the terminals of a reverse metering device are closed, the measured data is considered defective or an outlier. Of course, these two concepts are interrelated, and as stated above, by assuming a non-Gaussian distribution for noise measurements, defective data, or outliers, the likelihood of the improper performance of the nonrobust state estimation process is increased. Terefore, in addition to taking into account the non-Gaussian noise of the measurements, robust state estimation should be able to eliminate or reduce the adverse efect of defective data on the fnal output of the state estimation process.
Te robust state estimation (RSE) method used in this research has several advantages, which are as follows: (i) Consideration of non-Gaussian noise measurement (common positive feature with RMHE method in reference [6]) (ii) Ability to manage several types of undesirable data (iii) Te similarity of the implementation of this method with the WLS method (iv) Implementing a two-stage robust state estimator (v) Utilizing both traditional measurement sets and measurements from PMUs (vi) Possibility of using the two-stage approach presented in reference [11] for the robust dynamic state estimation method proposed in this research In the frst stage of the two-stage approach presented in reference [11], a generalized maximum probability estimator (GM estimator), which uses traditional SCADA measurements, is applied to make the state estimation process robust. Tis estimator includes the robust scale estimator, the projection statistics (PS) ( [48]), and the Huber model of the maximum probability-based estimator.
One way to limit non-Gaussian measurement noise is to extract a robust scale estimator. Te robust scale estimator, represented by the symbol s in the formulation of this paper, can be calculated in two ways and afects r i (mentioned in the objective function of the problem) and leads to the standardization of the measured noise. Te formulation of the two methods for calculating the coefcient s and how to enter it into the RSE problem based on the iteration algorithm is described in Section 3. It should be noted that by defning a robust scale estimator, it is no longer necessary to assume Gaussian noise for measurements in the iterative algorithm of the state estimation process; because the coefcient s indicates the unknown noise of the measurements [11].
Another way to limit non-Gaussian measurement noise in the generalized maximum probability estimator, in addition to extracting the robust scale estimator, is to use the Huber model of the maximum probability estimator, which is also referred to as the monotone maximum probability estimator [43,44].

Proposed Formulation
It is worth noting that the DSE model used in this paper is based on the fact that the state variables are similar to the measured quantities (of the magnitude and phase of the voltage). Tis fact is true in WAMS. Tis assumption is not true in traditional measurement systems. In such systems, a nonlinear measurement function communicates between state variables and measurement quantities. Terefore, the DSE model described does not apply to a measurement system consisting of purely traditional measurements.
It should also be noted that due to the inherent diferences between traditional measurements and PMU data (which are reported much faster and are equipped with accurate timestamps), it is not technically possible to add traditional measurements to WAMS. Te reverse is possible; by adding some PMU data to the SCADA system. In such cases, the described DSE model can be used. For this purpose, a traditional state estimate is performed in the presence of only traditional measurements. Its output, which includes the magnitude and phase of the voltages, is then considered as a set of phased data for the proposed WAMS-based DSE model [49].
On the other hand, the model considered for RSE in reference [11] has been carried out in two stages, the frst stage of which is considered in this research as the RSE method. Te frst stage mentioned in reference [11] uses SCADA measurements, and the second stage uses PMU measurements. However, the state estimation method mentioned in this reference is robust static state estimation and does not use predictions (which is the second set of DSE process inputs).
According to the above explanations, the model presented in this research for RDSE, while using the features of the DSE method mentioned in Section 2, uses the capabilities of the RSE method. Furthermore, this model presents a robust state estimation module that can manage the non-Gaussian noise of the measurement and reduces the adverse efect of outliers. Tis model is also able to eliminate the adverse efect of various disturbances in the power system on the predictions.
According to the explanation of Section 2, traditional measurements can also be included in the DSE process based on optimization. However, this matter requires a two-stage process. Te frst stage is based on iterative algorithms and uses only traditional measurements, and the second stage is based on optimization. In addition to the frst stage, the outputs, which are the magnitude and phase of voltages, also use measurements from PMUs.
Te choice of an iteration-based state estimation method for the frst stage is optional. Due to the shortcomings of the DSE method mentioned in this paper in reducing the adverse efect of outliers on the fnal output of the state estimation, in the proposed model of this research, the selected method of state estimation for the frst stage is the RSE method. Te results of numerical studies show that the proposed method, on the one hand, both in the case of normal network changes and in the case of sudden changes (in the event of disturbance), and on the other hand, despite the outlier in traditional measurements and taking into account the non-Gaussian noise, shows good performance.

Introduction of the Formulation of Dynamic State Estimation Method Based on Optimization.
Te method of calculating the indirect measurement values as, well as the predictions, are mentioned in the reference [3]. Te objective function and constraints of the optimized DSE method are defned as follows. It should be noted that the coefcients α i and β i in the Holt method are calculated through an optimization lateral problem in which a nonlinear optimization algorithm is used to minimize the mean square error index (MSE). Another method of calculating α i and β i coefcients is based on trial and error, which has not been used in this paper.

Te Objective Function of the Optimization Problem.
As mentioned in the initial description of this section about the introduction of the DSE method, in this method, the common state estimation problem based on the iteration algorithm becomes an optimization problem and thus, the amounts of power system state variables are estimated. First, the objective function of this optimization problem is defned in the following equation (1): In the aforementioned equation, x ijt , x it , and x it , respectively represent measurement j associated with the state variable i at period t, the estimation associated with state variable i at period t, and the prediction associated with state variable i at period t. Moreover, σ ij and σ ij are the standard deviation of measurement j associated with state variable i and the standard deviation of prediction associated with state variable i at period t, respectively. In the DSE method used in this paper, invalid predictions due to disturbances are discarded, which this discarding can be carried out in two ways; discarding by process (which can be solved by equating the J t derivative as the objective function with zero) and Simultaneous Discarding. In the second method, the objective function and constraints of the optimization problem are changed by applying a binary coefcient corresponding to invalid predictions. Terefore, the objective function of the problem will be in the form of the following equation (2): Te binary coefcient c it represents the sudden change in the state associated with the i-th state variable over the t-th period.

Constraints of the Optimization Problem.
Equations (3) and (4) are defned as follows and will be used in inequality (5): According to the issues raised in reference [3], the limitations of the binary coefcient c it are in the form of inequality (5).
M i is a positive constant number and as the fraction denominator is always assumed to be greater than the fraction numerator. When the fraction numerator is negative, the binary coefcient will be equal to zero, and when the fraction numerator is positive, the binary coefcient c it International Transactions on Electrical Energy Systems will be equal to 1. When the binary coefcient c it is equal to 1, the corresponding prediction data will be removed from the objective function.
Te problem that has been formed and addressed so far is nonlinear; because, frstly, in the objective function mentioned in equation (2), there is a multiplication of a binary variable and a component with power 2, and secondly, in relation to the auxiliary equation corresponding to the predictions (equation (3)), there is an absolute value of an expression. Terefore, the existing optimization problem is mixed integer nonlinear programming (MINLP). Conventional solvents do not have enough power to solve such problems efectively, even with a long execution time.

Formulation of Dynamic State Estimation Based on Mixed Integer
Programming. An efcient MIP solver can solve large-scale power system problems with numerous continuous, and integer variables in a tolerable amount of time. It also ensures a globally optimal solution or tolerable solution. MIP ofers highly fexible and accurate modeling capabilities that are highly desirable in practical problems. In the DSE method used in this research, the mixed-integer quadratic programming (MIQP) model is used. First, the equivalent variable R it is introduced according to equation (6).
In this way, the objective function mentioned in equation (2) is rewritten as equation (7).
Te above objective function is quadratic. To linearize the variable based on the variables and the inequalities mentioned in (8) and (9) are introduced [50].
If x it indicates the voltage magnitude, the value of r min i and r max i will be equal to 0 p.u. and 2 p.u., respectively, and if x it indicates the voltage angle, the value of r min i and r max i will be equal to − π and π radians, respectively. When the variable c it is equal to 1, the corresponding variable R it will be removed according to inequality (8) and the inequalities mentioned in (9) will be inefective. Conversely, when the variable c it is equal to 0, the corresponding variable R it is equal to r it according to inequality (9) and inequality (8) has no role in determining the variable R it . Tus, the nonlinear inequality (5) becomes the linear inequalities (8) and (9). Now, to linearize equations (3) and (4), new inequalities (10)-(13) are introduced.
Terefore, according to the process studied in this section, the state estimation problem is formulated in MIP. Te outputs of the aforementioned optimization problem are the estimated values for the system state variables (magnitude and phase of voltage of busbars).

Introduction of the Formulation of Robust State Estimation Method Based on Iteration
Algorithm. For a power system with n busbars, the relationship between the mdimensional measurement vector z and the state vector x (including the magnitude and phase of voltage of busbars) is described as an equation. To solve this state estimation problem, the common WLS method is used, which uses the Gauss-Newton iteration algorithm. At the same time, its performance is reduced in the presence of outliers or the absence of an error covariance matrix (due to non-Gaussian measurement noise). In this case, the potential features of the RSE method considered in this research (including the capabilities of the robust scale estimator and Huber model of maximum probability estimator) can be used [11].

Huber Model.
Non-Gaussian noise can follow wellknown distributions such as Laplace, Gaussian composition., or any other distribution. A person named Huber introduced a model of the maximum probability estimator (M estimator) in which the function is defned as equation (14) and is considered a convex nonlinear function. Tis function assigns a least absolute deviations (LAD) function (ρ(r) � |r|) to the outliers and an ordinary least squares (OLS) function (ρ(r) � r 2 ) to the other data according to equation (14).
In equation (14), c is a fxed number and r S i represents the standardized residual value and is defned as equation (15).
In equation (15), r i is the residual value (the diference between the function of the state variables to be estimated and the measured value of i) and is equal to r i � z i − h i (x). Also, s is the robust scale estimator. Coefcient ω is used in the objective function of the RSE problem and by using the method of projection statistics, it limits the efect of the leverage points. Te method of calculating ω i is explained as follows. 6 International Transactions on Electrical Energy Systems

Te Objective Function of the Problem.
To manage non-Gaussian noise measurement with optimal performance in the estimation process, along with making the state estimation robust to outliers, the objective function mentioned in equation (16) should be minimized.
Te weighting factor ω i , which limits the efect of the leverage points using the projection statistics method, is calculated through equation (17).
Te value of d i is calculated according to equation (18).
In the above equation, F represents the cumulative distribution function (CDF), V 1 is a positive numerical value between the values 0 and 1, and V 2 equals the sum of the nonzero values of line i of the Jacobin matrix of measurement (H).
In the above equation, Γ is the gamma distribution function. Assuming d i is fxed at d i � 1.5 [11], there will be a good performance for Gaussian noise. Of course, this coefcient needs to be adjusted to achieve the desired performance of state estimation for non-Gaussian noise. In this paper, the value of d i is calculated according to equations (18) and (19).

Calculating the Robust Scale Estimator.
Tere are two ways to calculate a robust scale estimator (s). Te frst method, which is more widely used, is the method called median absolute deviation (MAD). Te reason for the widespread use of this method is the simpler relationship and the possibility of faster calculation. Te defnition of s according to this method is based on equation (20).
Te correction factor b m is used to make a bias for Gaussian noise and must be adjusted for other symmetric non-Gaussian distributions to reduce the bias of the estimation. It is also important to note that s 2 tends to the variance of measurements (σ 2 ) when the number of samples to be estimated is infnite [51]. Although the MAD method is used in symmetric distributions, it has low-performance efciency under Gaussian noise. Tis can be improved by using the second method of calculating s in the following equation (21): In the above equation, f m is a correction factor and the lomed operator stands for the low median. Te outer lomed operator is the lower median (([m + 1]/2) th number among m numbers) and the inner lomed operator is the upper median ([m/2] +1-th number among m numbers). Te meaning of symbol [] in the interpretation of the ([m + 1]/2)-th number is the operator of absolute value. Te second method of calculating s ofers a better fnite sample efciency than the MAD method, even for completely asymmetric non-Gaussian distributions; while also maintaining the robustness of the estimation process [11].

Problem-Solving Algorithm.
To minimize the objective function mentioned in equation (16), the derivative J with respect to x must be set to zero.
In the above equation, a i is the i-th row of the Jacobin matrix (H � (zh/zx)) and ψ(r S i ) is equal to (zρ(r S i )/zr S i ). Te function ψ is a two-factor function and is obtained as follows (23): Multiplying and dividing on both sides of the above equation yields the matrix form mentioned in (24).
Te defnitions of Q and R correspond to equations (25) and (27).
Using the iterative-based algorithm called iteratively reweighted least squares (IRLS) and using the Taylor series frst-order expansion, correction of the state variable in k-iterations and the corresponding convergence limit for the robust state estimator considered in this paper will be according to equations (28) and (29), respectively.

International Transactions on Electrical Energy Systems
Starting the described fat voltage-based iteration algorithm (assuming the voltage magnitude is 1 and the voltage angle is 0 in the initial iteration) is not a good idea to initialize the proposed GM estimator; because in this case, most of the standardized residual values of the measurements are larger than the limit of the Huber function and are severely reduced in weight. Terefore, the most desirable solution is to perform WLS estimation for the frst iteration. Projection statistics calculation concerning fat voltage and the IRLS algorithm is then followed [11].

Introduction of the Formulation of a Two-Stage Robust Dynamic State Estimation Method Based on Iteration Algorithm and Optimization.
Te formulation of the frst and second stages of the RDSE method used in this paper can be achieved by applying some modifcations to the formulation of Sections 3.1 and 3.2. In the frst step, a generalized estimator based on the maximum probability should be used in terms of all its capabilities. Te output of the frst step is to determine the magnitude and phase of the voltage of all network busbars. Now, in the objective function, another part is added to the optimization-based dynamic state estimation method, which considers minimizing the diference between the outputs of the frst stage and the fnal estimated values for the state variables.
In the above equation, r RSE it is the diference between the output value of the frst stage corresponding to the variable of the i-th state in the time interval t and the corresponding value of the estimated i-th in the time interval t and σ RSE it is the amount of standard deviation of the output of the frst stage corresponding to the variable of the i-th state in the time interval t.
Te calculation method of r RSE it is similar to equations (3) and (4).
In the above equation, x RSE it represents the value of the output of the frst stage corresponding to the variable of the i-th state in the time interval t-th. Te above equation is also linearized using the new inequalities (32) and (33), and the absolute nonlinear operator is removed.
Te implementation of this two-stage process has the following advantages: (i) Ability to use both traditional measurements and measurements from PMUs (ii) Possibility of considering predictions along with measurements as two sets of DSE problem inputs (as one of the main features of DSE (iii) Ability to manage disturbances on the network and reduce their adverse efects (iv) Ability to make the state estimation process robust (v) Increase the redundancy of state estimation inputs by adding a third set of inputs to the problem (which are the same as the frst stage outputs), and consequently, improve the performance of the state estimator

Introducing the State Estimator Performance Evaluation
Index. To evaluate the performance of the state estimator, the mean square error index (MSE) is typically used in papers, which is defned according to (34).
In equation (34), NS is the number of state variables (or magnitude or phase of the voltage of the busbars), x i is the estimated value for the i-th state variable and x i is the actual value of the i-th state variable.
Te true values of state variables are the values obtained from the network load distribution and are available for IEEE standard networks. In a real network, the state estimation process itself is considered a prelude to optimal power fow (OPF). However, in this paper and similar studies, IEEE standard networks (where it is possible to obtain real values of state variables) are used to evaluate the performance of the proposed state estimator by comparing the values estimated by the state estimator with the expected values (obtained from power fow).
In this paper, according to the study of disturbances in the network at diferent times, the problem is investigated in several time intervals, and therefore, another performance evaluation index called the root mean square deviation (RMSD) index is used, which is defned in the following equation (35): Te RMSD value for the magnitude and phase of the voltages is calculated separately; because the scale of these two variables is diferent [2]. It is obvious that the lower the RMSD value, the better the performance of the state estimator is. In other words, the estimated values are more accurate. Conversely, the higher the RMSD value, the more unfavorable the performance of the state estimator in those circumstances is.

Numerical Studies
Te test networks considered in this paper are a 3-bus network, IEEE standard 14-bus network, and IEEE standard 57-bus network.
In examining the DSE method considered in this research (the second stage of the proposed RDSE method), the real values of the network state variables are obtained by performing power fow through the use of the MAT-POWER tool in MATLAB software. It is also assumed that direct measurement values (measurements derived from PMUs) are obtained by adding a random number in certain ranges to the corresponding real values. Te values of the predictions are also calculated by considering the explanations of the relevant section using the Holt method through the use of an optimization program in GAMS software and are entered into the main module of the program in MATLAB software. Te measurement error of the PMU based on the maximum allowable TVE is 1% [52]. For the uncertainty of transmission line parameters, the number 2% is assumed [53]. To apply diferent disturbances to the network, the DSE method has been implemented in several consecutive time intervals, and in some of these time intervals, diferent types of disturbances have been considered. Te number of time intervals is NT � 24.
In examining the RSE method considered in this research (the frst stage of the proposed RDSE method), the coefcient c (related to (3)) is assumed to be equal to 1.5, and the convergence limit of the iteration algorithm is assumed to be equal to 10 − 6 .
For each of the test networks, diferent types of scenarios are examined and analyzed to clarify the benefts of the ability to discard invalid prediction values in the event of network disturbances. Moreover, several disturbances have been assumed at diferent times to evaluate the performance of the proposed state estimator. In addition, diferent outliers have been hypothesized to investigate the efect of the robustness of the RDSE method.

3-Bus
Network. An overview of the 3-bus network under consideration is given in [27]. Te amount of base power is assumed as S base � 100MVA. Moreover, the production capacity of busbar 1 is equal to 0.8 p.u. + 0.25 p.u. . Busbar 2 is assumed to be the slack busbar and busbar 1 is the PMUequipped busbar. Te PMU in busbar 1 measures the phase value of bus voltage 1 and the current phase values of lines 1 and 2, and therefore the phase voltage values of busbars 2 and 3 are calculated indirectly. It is assumed that the bus load 3 changes at diferent time intervals. Busbar 3 load information and network line information are presented in Tables 1 and 2. Due to the direct measurement of the bus angle by the PMU, the reference bus angle is not considered. Te various assumptions for this network in the relevant studies are presented in Table 3.
When the generator of busbar 1 is disconnected from the network, the network load and network losses are supplied by busbar 2.

IEEE Standard 14-Bus
Network. An overview of the IEEE standard 14-bus network can be seen in [54]. Unnecessary mention of the specifcations of this network (including line information, increasing costs of power plant production, transformer tap changers, and busbar information) is avoided in this research. Tis information is available from various sources, including reference [55].
Buses 2, 6, 9, and 11 are assumed to be equipped with PMUs [54]. Terefore, the measurements of other busbars are obtained through the indirect measurement process. Table 4 presents the various assumptions for this network in the relevant studies.

Scenarios under Consideration.
To evaluate the performance of the state estimators mentioned in the previous sections, several scenarios are considered in this research. Te scenarios under consideration are as follows: (i) Scenario 1: implementation of SSE (by optimization method), by removing all prediction values from the analysis process and assuming no disturbance occurs (ii) Scenario 2: execute DSE (by optimization method), assuming no disturbance occurs (iii) Scenario 3: execute DSE (by optimization method), assuming a disturbance occurs, without discarding invalid predictions (iv) Scenario 4: execute DSE (by optimization method), assuming disturbance occurs, discarding invalid predictions (v) Scenario 5: execution of SSE using WLS (based on iteration algorithm) (vi) Scenario 6: RSE implementation (based on iteration algorithm) (vii) Scenario 7: implementation of two-stage DSE based on iteration and optimization algorithm (viii) Scenario 8: implementation of two-stage RDSE based on iteration and optimization algorithm Scenarios 1 to 6 refect the existing works, and scenarios 7 and 8 refect the proposed method of this paper. Diferent comparisons are conducted by presenting several tables and   fgures. It should be noted that scenarios 5 and 6 run at onetime interval, and the rest of the scenarios run at 24-time intervals. Terefore, the RMSD index is not used to evaluate the performance of scenarios 5 and 6.

Numerical Studies of Optimization-Based Dynamic State
Estimation Method. In this section, the performance of the state estimator under diferent conditions is compared. Specifcally, scenarios 1 to 4 are considered in this section. Te value of RMSD obtained in diferent cases is according to Table 6. Considering the results, the following points can be deduced: (i) Comparison of scenarios 1 and 2: calculating the predictions and adding them to the state estimation process, while increasing the redundancy of the input data to the state estimation process, improves (decreases) the RMSD. (ii) Comparison of scenarios 2 and 3: the occurrence of disturbance will increase the RMSD if the invalid predictions are not discarded. (iii) Comparison of scenarios 3 and 4: excluding invalid predictions in case of perturbation leads to an improvement (reduction) of RMSD. (iv) Comparison of scenarios 2 and 4: although disturbance increases the RMSD to some extent (undesirable), the RMSD can be improved by discarding invalid forecast data; however, it will not be the same as it was before the disturbance. Table 7 shows the values of binary coefcient c it for the 3-bus test network in scenario 4 at diferent times. Te numbers 0 and 1 in this table indicate that the corresponding prediction data are included and excluded, respectively. Table 8 presents the complete data of the busbar 3 phasor in the 3-bus network, including indirect measurement, calculated prediction, fnal estimated values, and actual values extracted from the power fow, with an accuracy of 4 decimal places in scenario 4.
Te following points are noteworthy in connection with Table 7: (i) Te value of the binary coefcient c it corresponding to the state variable V 2 (reference bus voltage magnitude) is equal to 0 at all times. (ii) Te value of the binary coefcient c it corresponding to the state variable V 3 at times 9, 20, and 23, corresponding to the type 2, 5, and 6 disturbances in Table 3, is equal to 1. Because busbar 3 is close to the site of the disturbance, and the corresponding amount of prediction is afected by the disturbance. (iii) Due to the type of calculation of the prediction (which is related to several previous time periods) and due to the wide changes in the operation of the power network in the assumed disturbances, the value of the binary coefcient c it corresponding to the state variable V 3 at other time intervals is equal to 1. (iv) Te degree of infuence of the amount of prediction data corresponding to the state variable δ 3 is greater than the perturbations in the network, compared to other state variables.        Tables 9 and 10 show the binary coefcient values equal to 1 for PMU-free buses in the IEEE standard 14-bus and 57bus networks in scenario 4 at diferent times. Similar to 3bus networks, PMU-equipped busbars contain more accurate measurements than PMU-free busbars, and the corresponding binary coefcient is often equal to 1. Terefore, Tables 9 and 10 focus on buses without PMU.
Concerning Table 9, the following can be mentioned: (i) At time t � 5, the bus load 13 increases abruptly. At this time, the prediction value corresponding to the variable of state variable of magnitude and angle of the voltage of busbars 9, 10, 11, 12, and 13 are afected. As expected, the busbars near the site of the disturbance were afected. (ii) At time t � 12, bus generator 2 is tripped out. Tis does not constitute a sudden change in the voltage magnitude state variables, and therefore the binary coefcient c it corresponding to none is not equal to 1 at t � 12. Tis is because the generator is close to the source of the disturbance (busbar 1 generator), and the busbar 3 capacitive compensator prevents a sharp voltage drop as the busbar 2 generator leaves the network. (iii) Te outage of the line between busbars 6 and 11 at time t � 18 leads to the discarding of the prediction data corresponding to the magnitude and phase of voltage variables of busbars 11, 12, and 13 at times t � 17 and t � 18, all of which somehow cause this disturbance. (iv) It is observed that the occurrence of a disturbance at a certain time also afects the prediction values of subsequent times, which is due to the fact that the values of previous time periods are taken into account in calculating the prediction value related to that specifc time.
Similar cases can be mentioned in connection with Table 10. It is only necessary to pay attention to the fact that the occurrence of disturbances in the 57-bus network has    less efect on the set of network state variables than the occurrence of disturbances in the 14-bus network due to the larger dimensions of the studied network.

Numerical Studies of Robust State Estimation Method
Based on Iteration Algorithm. In this section, scenarios 5 (WLS) and 6 (RSE) are examined and compared. For this purpose, we consider the IEEE standard 14-bus network and observe the role of RSE process robustness by applying diferent types of outliers. Te diagrams in Figure 1 show the output of the state estimation process for scenarios 5 and 6 without applying outliers.
Te diagrams in Figure 2 show the output of the state estimation process for scenarios 5 and 6 by applying an outlier. Te type of this piece of data is the change of the active power fow between busbars 2 and 3. When applying this outlier, it is assumed that the noise of the measurements follows the Gaussian distribution. Figure 1 and 2 diagrams illustrate the adverse efect of the presence of outliers on the fnal output of the WLS estimation process due to the high error in the estimated values and the role of the RSE process robustness in controlling this adverse efect.
Now it is assumed that the noise of the measurements is non-Gaussian, and the corresponding probability density  function is obtained from f(x) � (1/σ)e − ((x− μ) 2 /σ) . Te diagrams in Figure 3 show the output of the state estimation process for scenarios 5 and 6 under these conditions. Below each diagram, the value of the mean error value of the estimated values for the bus voltage magnitude is mentioned.
Te diagrams in Figure 3 show the increase in the error of the fnal estimated values of the state estimation process by implementing the WLS method (scenario 5). Due to the capabilities of the RSE method (scenario 6), including the use of a robust scale estimator and the Huber model of the maximum probability estimator, the adverse efect of non-Gaussian noise on the fnal output of the state estimation process is manageable.

Numerical Studies of Two-Stage Robust Dynamic State Estimation Method Based on Iteration Algorithm and
Optimization. In this section, by examining scenarios 7 and 8, the performance of the state estimator is compared. Table 11 presents the simulation results performed on the IEEE standard 14-bus network.

Sensitivity Analysis.
Te results show that the proposed RDSE method, while providing a good performance (low RMSD), can simultaneously manage the efect of network disturbances on predictions and the efect of outliers on the measurements. By examining Tables 6 and  11 and comparing scenarios 4 and 7, it is found that increasing the redundancy of the input data to the DSE process improves the performance of the state estimator. Tis improvement from the perspective of the RMSD m and RMSD p indices is equal to 21.59% and 18.37%, respectively. However, the presence of an outlier can still weaken the performance of the state estimator. Tis weakening from the perspective of the RMSD m and by assuming diferent outliers of changing the reactive power of one busbar, changing the active power of one busbar, and changing the active power of one line is equal to 63.77%, 40.58%, and 46.38%, respectively. Similarly, the weakening of the state estimation performance from the perspective of the RMSD p and by assuming the aforementioned diferent outliers is equal to 51.25%, 30%, and 40%, respectively. By considering the robustness feature in scenario 8 and comparing the obtained results from scenarios 7 and 8, the ability to control the adverse efect of outliers on the fnal output of the proposed two-stage    robust dynamic state estimation process is determined. From the perspective of the RMSD m index and considering simultaneously diferent disturbances, a non-Gaussian noise, and aforementioned diferent outliers, this improvement is equal to 25.66%, 26.8%, and 26.73%, respectively. Similarly, this improvement from the perspective of the RMSD p index and considering simultaneously the aforementioned conditions is equal to 24.79%, 25%, and 27.68, respectively. Figure 4 shows the output of the state estimation process for scenarios 7 and 8, assuming Gaussian and non-Gaussian noise measurements.
Examination of recent graphs confrms the benefcial efect of using the robustness of the state estimation feature in the proposed two-stage process; in this way, the adverse efect of non-Gaussian noise on the fnal outputs of the twostage process of state estimation is controlled, and the estimated values are closer to the true values.

Conclusions
In this paper, due to the importance of state estimation, this issue was addressed by considering its various aspects. In this regard, static and dynamic state estimations were compared. It was mentioned that dynamic state estimation has several capabilities, including the use of predictions along with measurements. Consequently, more redundancy of input data to the state estimation process and more accurate estimated outputs are provided. Due to the importance of dynamic state estimation in power systems nowadays, according to the capabilities of this type of state estimation, a method for dynamic state estimation in the presence of PMUs was investigated. Tis method was based on optimization and has a good performance against unwanted network disturbances. Furthermore, while the concept of robustness of the state estimation process was explained, a two-stage method was implemented to make the dynamic state estimation process robust against both the outliers and the non-Gaussian noise of measurements. Finally, a two-stage method named two-stage robust dynamic state estimation was proposed and thoroughly analyzed. Te efectiveness of this proposed method was verifed by assuming diferent kinds of probable disturbances such as load fuctuations and generator outages and diferent types of noise (Gaussian and non-Gaussian), along with the undesirable injection of outliers into three sample networks. Tis method increases the accuracy and enhances the redundancy of the state estimation process. By defning the root mean square deviation (RMSD) index, which in essence and based on its defnition, encompasses the defnition of two other usual indices; standard deviation and mean square error, diferent scenarios were introduced to compare the fnal output of the dynamic state estimator which is both robust and capable of disturbance rejection with other types of state estimators. It was observed that, for example, for the IEEE standard 14-bus network, the performance of the state estimation process from the perspective of the RMSD m index was improved to the extent of 4.55%, 19.32%, and 15.91%. Tese three fgures are for consideration of simultaneous diferent disturbances, a non-Gaussian noise, and diferent outliers of changing the reactive power of one busbar, changing the active power of one busbar, and changing the active power of one line, respectively. Similarly, the performance of the state estimation process from the perspective of the RMSD p index was improved to the extent of 7.14%, 20.41%, and 17.35%, considering simultaneously the aforementioned conditions.

Data Availability
Te data used to support the fndings of this study are included within the article.