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

Authors: Vahid Hakimian; Hossein Seifi
Journal: International Transactions on Electrical Energy Systems
 (2023)

---

## Abstract

State estimation provides the best estimate of the system states using a set of measurements that is divided into two main classes; static and dynamic types. The state estimator should provide an acceptable response in the presence of system changes and/or disturbances. To achieve this aim, considering the high rate of data changes and the high volume of exchanging information, the phasor measurement units could play an outstanding role. Furthermore, due to the inevitable errors of the measurements, and to provide an authentic estimate of the system states using the measured data, a state estimator should be able to eliminate or decrease the undesirable effect of bad data or outliers on the final output of the process. This type of state estimator is called a robust type. In this paper, a dynamic state estimator for a system with the synchronized phasor measurement units is utilized, which is capable of decreasing the effect of disturbances. Simultaneously, the robustness of the state estimation process in the presence of outliers is considered, which enables the state estimation process to reduce the negative effect of outliers on its final output. Therefore, the implemented process is called a robust dynamic state estimation method. By defining different scenarios and utilizing numerical indices, simulation results are thoroughly analyzed based on an explanatory three-bus network and the IEEE standard 14-bus and 57-bus networks using the MATPOWER tool.


---

## Body

## 1. Introduction

The 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 reflects the behavior of that system. The power system can be in any of the three states, normal state, emergency state, and restorative state. The 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]. There 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. The 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. The measurements must be synchronized chronologically to obtain an accurate picture of the power system [3].The 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. The 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. The information received from the RTU is processed by a state estimator [4]. The 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 fiber 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 intensification of the influence of phasor measurement units (PMUs). These capabilities include more accurate measurements from PMUs, higher sampling rates, temporal synchronization of data, and linearity between measurement values and state variables. Thus, the effect 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 field 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 different approaches to the subject of state estimation of PMUs, which fall into several main categories of modification 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–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–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. This 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. The 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 different kinds of network sudden changes (disturbances). The 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 filter, this method has some advantages, such as the lack of dependence on system linearity and Gaussian noise, resistance to falsified 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 first stage of which uses a generalized maximum probability estimator (GM estimator). This GM estimator uses traditional SCADA measurements. In the second stage, the results of the first 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–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 different control modes, and [35] estimates the state in the presence of FACTS and devices of unified power flow control (UPFC) type. [36–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. This 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, definitions, 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–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 effective 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 effect of such defective or bad data on the final output of the process which is called the robust state estimator. The lack of enough related research at the transmission level of power systems considering the simultaneous significant effects 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 effect of disturbances on the network, the robustness of the state estimation process is also considered against the defective data affecting the existing measurements. Moreover, unlike the conventional repetitive algorithms and methods, an optimization-based 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. This 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 confidence bound. By implementing an explanatory three-bus network and the IEEE standard 14-bus and 57-bus networks using the MATPOWER tool, simulation results are thoroughly analyzed.The rest of this paper is presented as follows. Section2 provides an overview of the DSE and RSE. The proposed RDSE model and formulation are presented in Section 3. Section 4 discusses numerical examples, and the final points are stated in Section 5.


## 2. 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. Thus, 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 filter are used to estimate the system states. The accuracy of the DSE depends on the sampling rate of the measurements. The 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. The introduction of the PMU with a high sampling rate between 30–60 Hz made it possible to implement the Kalman filter algorithm to implement DSE [4].The power system follows a quasistatic behavior under normal operating conditions; the daily profile is subject to small random fluctuations. Under these conditions, the power system state changes slowly, and the DSE performance is acceptable. However, during various disturbances (such as load fluctuations, generator outages, and network switchings), the state variables corresponding to the busbars that are close to the area affected 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 effect of the invalid predictions on the output values of the estimation state will reduce the performance of the dynamic estimation process.The 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. This 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 effect 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 specific 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 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]. The Holt method used in this section is a conventional representation of the Holt method used in several kinds of research.


### 2.2. 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 affected by deviations from the true values of several measurements, then the corresponding estimator is considered robust. The 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 effect 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]. The 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 effect 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 different concepts; i.e., the first case shows the difference 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. Therefore, 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 effect of defective data on the final output of the state estimation process.The 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)
The 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 researchIn the first 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. This 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. The robust scale estimator, represented by the symbols in the formulation of this paper, can be calculated in two ways and affects ri (mentioned in the objective function of the problem) and leads to the standardization of the measured noise. The formulation of the two methods for calculating the coefficient 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 defining 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 coefficient 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].


## 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. Thus, 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 filter are used to estimate the system states. The accuracy of the DSE depends on the sampling rate of the measurements. The 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. The introduction of the PMU with a high sampling rate between 30–60 Hz made it possible to implement the Kalman filter algorithm to implement DSE [4].The power system follows a quasistatic behavior under normal operating conditions; the daily profile is subject to small random fluctuations. Under these conditions, the power system state changes slowly, and the DSE performance is acceptable. However, during various disturbances (such as load fluctuations, generator outages, and network switchings), the state variables corresponding to the busbars that are close to the area affected 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 effect of the invalid predictions on the output values of the estimation state will reduce the performance of the dynamic estimation process.The 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. This 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 effect 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 specific 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 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]. The Holt method used in this section is a conventional representation of the Holt method used in several kinds of research.


## 2.2. 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 affected by deviations from the true values of several measurements, then the corresponding estimator is considered robust. The 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 effect 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]. The 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 effect 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 different concepts; i.e., the first case shows the difference 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. Therefore, 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 effect of defective data on the final output of the state estimation process.The 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)
The 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 researchIn the first 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. This 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. The robust scale estimator, represented by the symbols in the formulation of this paper, can be calculated in two ways and affects ri (mentioned in the objective function of the problem) and leads to the standardization of the measured noise. The formulation of the two methods for calculating the coefficient 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 defining 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 coefficient 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].


## 3. 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). This fact is true in WAMS. This assumption is not true in traditional measurement systems. In such systems, a nonlinear measurement function communicates between state variables and measurement quantities. Therefore, 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 differences 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. The 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 first stage of which is considered in this research as the RSE method. The first 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 Section2, 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 effect of outliers. This model is also able to eliminate the adverse effect of various disturbances in the power system on the predictions.According to the explanation of Section2, traditional measurements can also be included in the DSE process based on optimization. However, this matter requires a two-stage process. The first stage is based on iterative algorithms and uses only traditional measurements, and the second stage is based on optimization. In addition to the first stage, the outputs, which are the magnitude and phase of voltages, also use measurements from PMUs.The choice of an iteration-based state estimation method for the first stage is optional. Due to the shortcomings of the DSE method mentioned in this paper in reducing the adverse effect of outliers on the final output of the state estimation, in the proposed model of this research, the selected method of state estimation for the first stage is the RSE method. The 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.


### 3.1. Introduction of the Formulation of Dynamic State Estimation Method Based on Optimization

The method of calculating the indirect measurement values as, well as the predictions, are mentioned in the reference [3]. The objective function and constraints of the optimized DSE method are defined as follows. It should be noted that the coefficients α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 coefficients is based on trial and error, which has not been used in this paper.


#### 3.1.1. The 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 defined in the following equation (1):(1)minJt=∑i=1NS∑j=1NMix¯ijt−x^it2σ¯ij2+∑i=1NSx˜it−x^it2σ˜ij2.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 Jt 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 coefficient corresponding to invalid predictions. Therefore, the objective function of the problem will be in the form of the following equation (2):(2)minJt=∑i=1NS∑j=1NMix¯ijt−x^it2σ¯ij2+∑i=1NS1−γitx˜it−x^it2σ˜ij2.The binary coefficientγit represents the sudden change in the state associated with the i-th state variable over the t-th period.


#### 3.1.2. Constraints of the Optimization Problem

Equations (3) and (4) are defined as follows and will be used in inequality (5):(3)r˜it=x˜it−x^it,(4)r¯ijt=x¯ijt−x^it.According to the issues raised in reference [3], the limitations of the binary coefficient γit are in the form of inequality (5).(5)r˜it−3σ^itMi≤γit≤r˜it−3σ^itMi+1.Mi 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 coefficient will be equal to zero, and when the fraction numerator is positive, the binary coefficient γit will be equal to 1. When the binary coefficient γit is equal to 1, the corresponding prediction data will be removed from the objective function.The problem that has been formed and addressed so far is nonlinear; because, firstly, 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. Therefore, the existing optimization problem is mixed integer nonlinear programming (MINLP). Conventional solvents do not have enough power to solve such problems effectively, even with a long execution time.


#### 3.1.3. Formulation of Dynamic State Estimation Based on Mixed Integer Programming

An efficient 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 offers highly flexible 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 variableR˜it is introduced according to equation (6).(6)R˜it2=1−γitr˜it2.In this way, the objective function mentioned in equation (2) is rewritten as equation (7).(7)minJt=∑i=1NS∑j=1NMir¯ijt2σ¯ij2+∑i=1NSR˜it2σ˜it2.The above objective function is quadratic. To linearize the variable based on the variables and the inequalities mentioned in (8) and (9) are introduced [50].(8)rimin1−γit≤R˜it≤rimax1−γit,(9)r˜it−rimaxγit≤R˜it≤r˜it−riminγit.Ifx^it indicates the voltage magnitude, the value of rimin and rimax will be equal to 0p.u. and 2p.u., respectively, and if x^it indicates the voltage angle, the value of rimin and rimax will be equal to −π and π radians, respectively. When the variable γ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 ineffective. Conversely, when the variable γ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. Thus, the nonlinear inequality (5) becomes the linear inequalities (8) and (9).Now, to linearize equations (3) and (4), new inequalities (10)–(13) are introduced.(10)r¯ijt≥x¯ijt−x^it,(11)r¯ijt≥x^it−x¯ijt,(12)r˜it≥x˜it−x^it,(13)r˜it≥x^it−x˜it.Therefore, according to the process studied in this section, the state estimation problem is formulated in MIP. The outputs of the aforementioned optimization problem are the estimated values for the system state variables (magnitude and phase of voltage of busbars).


### 3.2. Introduction of the Formulation of Robust State Estimation Method Based on Iteration Algorithm

For a power system withn busbars, the relationship between the m-dimensional 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].


#### 3.2.1. Huber Model

Non-Gaussian noise can follow well-known 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 defined as equation (14) and is considered a convex nonlinear function. This function assigns a least absolute deviations (LAD) function ρr=r to the outliers and an ordinary least squares (OLS) function ρr=r2 to the other data according to equation (14).(14)ρrSi=12rSi2,forrSi<c,crSi−c22,elsewhere.In equation (14), c is a fixed number and rSi represents the standardized residual value and is defined as equation (15).(15)rSi=risωi.In equation (15), ri is the residual value (the difference between the function of the state variables to be estimated and the measured value of i) and is equal to ri=zi−hix. Also, s is the robust scale estimator. Coefficient ω is used in the objective function of the RSE problem and by using the method of projection statistics, it limits the effect of the leverage points. The method of calculating ωi is explained as follows.


#### 3.2.2. The 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.(16)J=∑i=1mωi2ρrSi.The weighting factorωi, which limits the effect of the leverage points using the projection statistics method, is calculated through equation (17).(17)ωi=min1,diPSi2.The value ofdi is calculated according to equation (18).(18)di=χ2−1V1,V2=F−1V1V2=di:FdiV2=V1.In the above equation,F represents the cumulative distribution function (CDF), V1 is a positive numerical value between the values 0 and 1, and V2 equals the sum of the nonzero values of line i of the Jacobin matrix of measurement (H).Concerning equations (18), and (19) is established.(19)V1=FdiV2=∫0xtV2−2/2e−t/22V2/2ΓV2/2dt.In the above equation, Г is the gamma distribution function. Assumingdi is fixed at di = 1.5 [11], there will be a good performance for Gaussian noise. Of course, this coefficient needs to be adjusted to achieve the desired performance of state estimation for non-Gaussian noise. In this paper, the value of di is calculated according to equations (18) and (19).


#### 3.2.3. Calculating the Robust Scale Estimator

There are two ways to calculate a robust scale estimator (s). The first method, which is more widely used, is the method called median absolute deviation (MAD). The reason for the widespread use of this method is the simpler relationship and the possibility of faster calculation. The definition of s according to this method is based on equation (20).(20)s=1.4826∙bm∙medianiri−medianjrj.The correction factorbm 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 s2 tends to the variance of measurements (σ2) when the number of samples to be estimated is infinite [51]. Although the MAD method is used in symmetric distributions, it has low-performance efficiency under Gaussian noise. This can be improved by using the second method of calculating s in the following equation (21):(21)s=1.1926⋅fm⋅lomedlomedi=1,⋯,mj≠iri−rj.In the above equation,fm is a correction factor and the lomed operator stands for the low median. The 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). The meaning of symbol [] in the interpretation of the m+1/2-th number is the operator of absolute value. The second method of calculating s offers a better finite sample efficiency than the MAD method, even for completely asymmetric non-Gaussian distributions; while also maintaining the robustness of the estimation process [11].


#### 3.2.4. Problem-Solving Algorithm

To minimize the objective function mentioned in equation (16), the derivative J with respect to x must be set to zero.(22)∂J∂x=∑i=1m−ωiaisi2ψrSi=0.In the above equation,ai is the i-th row of the Jacobin matrix (H=∂h/∂x) and ψrSi is equal to ∂ρrSi/∂rSi. The function ψ is a two-factor function and is obtained as follows (23):(23)ψrSi=rSi;rSi<c,ψrSi=csignrSi;rSi≥c.Multiplying and dividing on both sides of the above equation yields the matrix form mentioned in (24).(24)HTQR^−1z−hx=0.The definitions ofQ and R^ correspond to equations (25) and (27).(25)Q=diagqrSi,(26)qrSi=ψrSirSi,(27)R^=diagsi2.Using the iterative-based algorithm called iteratively reweighted least squares (IRLS) and using the Taylor series first-order expansion, correction of the state variable ink-iterations and the corresponding convergence limit for the robust state estimator considered in this paper will be according to equations (28) and (29), respectively.(28)HTQkR^−1H∆xk=HTQkR^−1z−hxk,(29)∆xk∞=xk+1−xk∞≤10−6.Starting the described flat 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. Therefore, the most desirable solution is to perform WLS estimation for the first iteration. Projection statistics calculation concerning flat voltage and the IRLS algorithm is then followed [11].


### 3.3. Introduction of the Formulation of a Two-Stage Robust Dynamic State Estimation Method Based on Iteration Algorithm and Optimization

The formulation of the first and second stages of the RDSE method used in this paper can be achieved by applying some modifications to the formulation of Sections 3.1 and 3.2. In the first step, a generalized estimator based on the maximum probability should be used in terms of all its capabilities. The output of the first 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 difference between the outputs of the first stage and the final estimated values for the state variables.(30)MinJt=∑i=1NSrRSEit2σRSEit2+∑i=1NS∑j=1NMir¯ijt2σ¯ij2+∑i=1NSR˜it2σ˜it2.In the above equation,rRSEit is the difference between the output value of the first 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 σRSEit is the amount of standard deviation of the output of the first stage corresponding to the variable of the i-th state in the time interval t.The calculation method ofrRSEit is similar to equations (3) and (4).(31)rRSEit=xRSEit−x^it.In the above equation,xRSEit represents the value of the output of the first stage corresponding to the variable of the i-th state in the time interval t-th. The above equation is also linearized using the new inequalities (32) and (33), and the absolute nonlinear operator is removed.(32)rRSEit≥xRSEit−x^it,(33)rRSEit≥x^it−xRSEit.The 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 effects(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 first stage outputs), and consequently, improve the performance of the state estimator


### 3.4. 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 defined according to (34).(34)MSE=1NS∑i=1NSx^i−xi2.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 xi is the actual value of the i-th state variable.The 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 flow (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 flow).In this paper, according to the study of disturbances in the network at different 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 defined in the following equation (35):(35)RMSD=1NT∑t=1NT1NS∑i=1NSx^it−xit2.The RMSD value for the magnitude and phase of the voltages is calculated separately; because the scale of these two variables is different [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.


## 3.1. Introduction of the Formulation of Dynamic State Estimation Method Based on Optimization

The method of calculating the indirect measurement values as, well as the predictions, are mentioned in the reference [3]. The objective function and constraints of the optimized DSE method are defined as follows. It should be noted that the coefficients α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 coefficients is based on trial and error, which has not been used in this paper.


### 3.1.1. The 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 defined in the following equation (1):(1)minJt=∑i=1NS∑j=1NMix¯ijt−x^it2σ¯ij2+∑i=1NSx˜it−x^it2σ˜ij2.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 Jt 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 coefficient corresponding to invalid predictions. Therefore, the objective function of the problem will be in the form of the following equation (2):(2)minJt=∑i=1NS∑j=1NMix¯ijt−x^it2σ¯ij2+∑i=1NS1−γitx˜it−x^it2σ˜ij2.The binary coefficientγit represents the sudden change in the state associated with the i-th state variable over the t-th period.


### 3.1.2. Constraints of the Optimization Problem

Equations (3) and (4) are defined as follows and will be used in inequality (5):(3)r˜it=x˜it−x^it,(4)r¯ijt=x¯ijt−x^it.According to the issues raised in reference [3], the limitations of the binary coefficient γit are in the form of inequality (5).(5)r˜it−3σ^itMi≤γit≤r˜it−3σ^itMi+1.Mi 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 coefficient will be equal to zero, and when the fraction numerator is positive, the binary coefficient γit will be equal to 1. When the binary coefficient γit is equal to 1, the corresponding prediction data will be removed from the objective function.The problem that has been formed and addressed so far is nonlinear; because, firstly, 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. Therefore, the existing optimization problem is mixed integer nonlinear programming (MINLP). Conventional solvents do not have enough power to solve such problems effectively, even with a long execution time.


### 3.1.3. Formulation of Dynamic State Estimation Based on Mixed Integer Programming

An efficient 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 offers highly flexible 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 variableR˜it is introduced according to equation (6).(6)R˜it2=1−γitr˜it2.In this way, the objective function mentioned in equation (2) is rewritten as equation (7).(7)minJt=∑i=1NS∑j=1NMir¯ijt2σ¯ij2+∑i=1NSR˜it2σ˜it2.The above objective function is quadratic. To linearize the variable based on the variables and the inequalities mentioned in (8) and (9) are introduced [50].(8)rimin1−γit≤R˜it≤rimax1−γit,(9)r˜it−rimaxγit≤R˜it≤r˜it−riminγit.Ifx^it indicates the voltage magnitude, the value of rimin and rimax will be equal to 0p.u. and 2p.u., respectively, and if x^it indicates the voltage angle, the value of rimin and rimax will be equal to −π and π radians, respectively. When the variable γ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 ineffective. Conversely, when the variable γ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. Thus, the nonlinear inequality (5) becomes the linear inequalities (8) and (9).Now, to linearize equations (3) and (4), new inequalities (10)–(13) are introduced.(10)r¯ijt≥x¯ijt−x^it,(11)r¯ijt≥x^it−x¯ijt,(12)r˜it≥x˜it−x^it,(13)r˜it≥x^it−x˜it.Therefore, according to the process studied in this section, the state estimation problem is formulated in MIP. The outputs of the aforementioned optimization problem are the estimated values for the system state variables (magnitude and phase of voltage of busbars).


## 3.1.1. The 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 defined in the following equation (1):(1)minJt=∑i=1NS∑j=1NMix¯ijt−x^it2σ¯ij2+∑i=1NSx˜it−x^it2σ˜ij2.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 Jt 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 coefficient corresponding to invalid predictions. Therefore, the objective function of the problem will be in the form of the following equation (2):(2)minJt=∑i=1NS∑j=1NMix¯ijt−x^it2σ¯ij2+∑i=1NS1−γitx˜it−x^it2σ˜ij2.The binary coefficientγit represents the sudden change in the state associated with the i-th state variable over the t-th period.


## 3.1.2. Constraints of the Optimization Problem

Equations (3) and (4) are defined as follows and will be used in inequality (5):(3)r˜it=x˜it−x^it,(4)r¯ijt=x¯ijt−x^it.According to the issues raised in reference [3], the limitations of the binary coefficient γit are in the form of inequality (5).(5)r˜it−3σ^itMi≤γit≤r˜it−3σ^itMi+1.Mi 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 coefficient will be equal to zero, and when the fraction numerator is positive, the binary coefficient γit will be equal to 1. When the binary coefficient γit is equal to 1, the corresponding prediction data will be removed from the objective function.The problem that has been formed and addressed so far is nonlinear; because, firstly, 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. Therefore, the existing optimization problem is mixed integer nonlinear programming (MINLP). Conventional solvents do not have enough power to solve such problems effectively, even with a long execution time.


## 3.1.3. Formulation of Dynamic State Estimation Based on Mixed Integer Programming

An efficient 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 offers highly flexible 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 variableR˜it is introduced according to equation (6).(6)R˜it2=1−γitr˜it2.In this way, the objective function mentioned in equation (2) is rewritten as equation (7).(7)minJt=∑i=1NS∑j=1NMir¯ijt2σ¯ij2+∑i=1NSR˜it2σ˜it2.The above objective function is quadratic. To linearize the variable based on the variables and the inequalities mentioned in (8) and (9) are introduced [50].(8)rimin1−γit≤R˜it≤rimax1−γit,(9)r˜it−rimaxγit≤R˜it≤r˜it−riminγit.Ifx^it indicates the voltage magnitude, the value of rimin and rimax will be equal to 0p.u. and 2p.u., respectively, and if x^it indicates the voltage angle, the value of rimin and rimax will be equal to −π and π radians, respectively. When the variable γ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 ineffective. Conversely, when the variable γ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. Thus, the nonlinear inequality (5) becomes the linear inequalities (8) and (9).Now, to linearize equations (3) and (4), new inequalities (10)–(13) are introduced.(10)r¯ijt≥x¯ijt−x^it,(11)r¯ijt≥x^it−x¯ijt,(12)r˜it≥x˜it−x^it,(13)r˜it≥x^it−x˜it.Therefore, according to the process studied in this section, the state estimation problem is formulated in MIP. The outputs of the aforementioned optimization problem are the estimated values for the system state variables (magnitude and phase of voltage of busbars).


## 3.2. Introduction of the Formulation of Robust State Estimation Method Based on Iteration Algorithm

For a power system withn busbars, the relationship between the m-dimensional 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].


### 3.2.1. Huber Model

Non-Gaussian noise can follow well-known 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 defined as equation (14) and is considered a convex nonlinear function. This function assigns a least absolute deviations (LAD) function ρr=r to the outliers and an ordinary least squares (OLS) function ρr=r2 to the other data according to equation (14).(14)ρrSi=12rSi2,forrSi<c,crSi−c22,elsewhere.In equation (14), c is a fixed number and rSi represents the standardized residual value and is defined as equation (15).(15)rSi=risωi.In equation (15), ri is the residual value (the difference between the function of the state variables to be estimated and the measured value of i) and is equal to ri=zi−hix. Also, s is the robust scale estimator. Coefficient ω is used in the objective function of the RSE problem and by using the method of projection statistics, it limits the effect of the leverage points. The method of calculating ωi is explained as follows.


### 3.2.2. The 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.(16)J=∑i=1mωi2ρrSi.The weighting factorωi, which limits the effect of the leverage points using the projection statistics method, is calculated through equation (17).(17)ωi=min1,diPSi2.The value ofdi is calculated according to equation (18).(18)di=χ2−1V1,V2=F−1V1V2=di:FdiV2=V1.In the above equation,F represents the cumulative distribution function (CDF), V1 is a positive numerical value between the values 0 and 1, and V2 equals the sum of the nonzero values of line i of the Jacobin matrix of measurement (H).Concerning equations (18), and (19) is established.(19)V1=FdiV2=∫0xtV2−2/2e−t/22V2/2ΓV2/2dt.In the above equation, Г is the gamma distribution function. Assumingdi is fixed at di = 1.5 [11], there will be a good performance for Gaussian noise. Of course, this coefficient needs to be adjusted to achieve the desired performance of state estimation for non-Gaussian noise. In this paper, the value of di is calculated according to equations (18) and (19).


### 3.2.3. Calculating the Robust Scale Estimator

There are two ways to calculate a robust scale estimator (s). The first method, which is more widely used, is the method called median absolute deviation (MAD). The reason for the widespread use of this method is the simpler relationship and the possibility of faster calculation. The definition of s according to this method is based on equation (20).(20)s=1.4826∙bm∙medianiri−medianjrj.The correction factorbm 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 s2 tends to the variance of measurements (σ2) when the number of samples to be estimated is infinite [51]. Although the MAD method is used in symmetric distributions, it has low-performance efficiency under Gaussian noise. This can be improved by using the second method of calculating s in the following equation (21):(21)s=1.1926⋅fm⋅lomedlomedi=1,⋯,mj≠iri−rj.In the above equation,fm is a correction factor and the lomed operator stands for the low median. The 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). The meaning of symbol [] in the interpretation of the m+1/2-th number is the operator of absolute value. The second method of calculating s offers a better finite sample efficiency than the MAD method, even for completely asymmetric non-Gaussian distributions; while also maintaining the robustness of the estimation process [11].


### 3.2.4. Problem-Solving Algorithm

To minimize the objective function mentioned in equation (16), the derivative J with respect to x must be set to zero.(22)∂J∂x=∑i=1m−ωiaisi2ψrSi=0.In the above equation,ai is the i-th row of the Jacobin matrix (H=∂h/∂x) and ψrSi is equal to ∂ρrSi/∂rSi. The function ψ is a two-factor function and is obtained as follows (23):(23)ψrSi=rSi;rSi<c,ψrSi=csignrSi;rSi≥c.Multiplying and dividing on both sides of the above equation yields the matrix form mentioned in (24).(24)HTQR^−1z−hx=0.The definitions ofQ and R^ correspond to equations (25) and (27).(25)Q=diagqrSi,(26)qrSi=ψrSirSi,(27)R^=diagsi2.Using the iterative-based algorithm called iteratively reweighted least squares (IRLS) and using the Taylor series first-order expansion, correction of the state variable ink-iterations and the corresponding convergence limit for the robust state estimator considered in this paper will be according to equations (28) and (29), respectively.(28)HTQkR^−1H∆xk=HTQkR^−1z−hxk,(29)∆xk∞=xk+1−xk∞≤10−6.Starting the described flat 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. Therefore, the most desirable solution is to perform WLS estimation for the first iteration. Projection statistics calculation concerning flat voltage and the IRLS algorithm is then followed [11].


## 3.2.1. Huber Model

Non-Gaussian noise can follow well-known 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 defined as equation (14) and is considered a convex nonlinear function. This function assigns a least absolute deviations (LAD) function ρr=r to the outliers and an ordinary least squares (OLS) function ρr=r2 to the other data according to equation (14).(14)ρrSi=12rSi2,forrSi<c,crSi−c22,elsewhere.In equation (14), c is a fixed number and rSi represents the standardized residual value and is defined as equation (15).(15)rSi=risωi.In equation (15), ri is the residual value (the difference between the function of the state variables to be estimated and the measured value of i) and is equal to ri=zi−hix. Also, s is the robust scale estimator. Coefficient ω is used in the objective function of the RSE problem and by using the method of projection statistics, it limits the effect of the leverage points. The method of calculating ωi is explained as follows.


## 3.2.2. The 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.(16)J=∑i=1mωi2ρrSi.The weighting factorωi, which limits the effect of the leverage points using the projection statistics method, is calculated through equation (17).(17)ωi=min1,diPSi2.The value ofdi is calculated according to equation (18).(18)di=χ2−1V1,V2=F−1V1V2=di:FdiV2=V1.In the above equation,F represents the cumulative distribution function (CDF), V1 is a positive numerical value between the values 0 and 1, and V2 equals the sum of the nonzero values of line i of the Jacobin matrix of measurement (H).Concerning equations (18), and (19) is established.(19)V1=FdiV2=∫0xtV2−2/2e−t/22V2/2ΓV2/2dt.In the above equation, Г is the gamma distribution function. Assumingdi is fixed at di = 1.5 [11], there will be a good performance for Gaussian noise. Of course, this coefficient needs to be adjusted to achieve the desired performance of state estimation for non-Gaussian noise. In this paper, the value of di is calculated according to equations (18) and (19).


## 3.2.3. Calculating the Robust Scale Estimator

There are two ways to calculate a robust scale estimator (s). The first method, which is more widely used, is the method called median absolute deviation (MAD). The reason for the widespread use of this method is the simpler relationship and the possibility of faster calculation. The definition of s according to this method is based on equation (20).(20)s=1.4826∙bm∙medianiri−medianjrj.The correction factorbm 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 s2 tends to the variance of measurements (σ2) when the number of samples to be estimated is infinite [51]. Although the MAD method is used in symmetric distributions, it has low-performance efficiency under Gaussian noise. This can be improved by using the second method of calculating s in the following equation (21):(21)s=1.1926⋅fm⋅lomedlomedi=1,⋯,mj≠iri−rj.In the above equation,fm is a correction factor and the lomed operator stands for the low median. The 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). The meaning of symbol [] in the interpretation of the m+1/2-th number is the operator of absolute value. The second method of calculating s offers a better finite sample efficiency than the MAD method, even for completely asymmetric non-Gaussian distributions; while also maintaining the robustness of the estimation process [11].


## 3.2.4. Problem-Solving Algorithm

To minimize the objective function mentioned in equation (16), the derivative J with respect to x must be set to zero.(22)∂J∂x=∑i=1m−ωiaisi2ψrSi=0.In the above equation,ai is the i-th row of the Jacobin matrix (H=∂h/∂x) and ψrSi is equal to ∂ρrSi/∂rSi. The function ψ is a two-factor function and is obtained as follows (23):(23)ψrSi=rSi;rSi<c,ψrSi=csignrSi;rSi≥c.Multiplying and dividing on both sides of the above equation yields the matrix form mentioned in (24).(24)HTQR^−1z−hx=0.The definitions ofQ and R^ correspond to equations (25) and (27).(25)Q=diagqrSi,(26)qrSi=ψrSirSi,(27)R^=diagsi2.Using the iterative-based algorithm called iteratively reweighted least squares (IRLS) and using the Taylor series first-order expansion, correction of the state variable ink-iterations and the corresponding convergence limit for the robust state estimator considered in this paper will be according to equations (28) and (29), respectively.(28)HTQkR^−1H∆xk=HTQkR^−1z−hxk,(29)∆xk∞=xk+1−xk∞≤10−6.Starting the described flat 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. Therefore, the most desirable solution is to perform WLS estimation for the first iteration. Projection statistics calculation concerning flat voltage and the IRLS algorithm is then followed [11].


## 3.3. Introduction of the Formulation of a Two-Stage Robust Dynamic State Estimation Method Based on Iteration Algorithm and Optimization

The formulation of the first and second stages of the RDSE method used in this paper can be achieved by applying some modifications to the formulation of Sections 3.1 and 3.2. In the first step, a generalized estimator based on the maximum probability should be used in terms of all its capabilities. The output of the first 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 difference between the outputs of the first stage and the final estimated values for the state variables.(30)MinJt=∑i=1NSrRSEit2σRSEit2+∑i=1NS∑j=1NMir¯ijt2σ¯ij2+∑i=1NSR˜it2σ˜it2.In the above equation,rRSEit is the difference between the output value of the first 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 σRSEit is the amount of standard deviation of the output of the first stage corresponding to the variable of the i-th state in the time interval t.The calculation method ofrRSEit is similar to equations (3) and (4).(31)rRSEit=xRSEit−x^it.In the above equation,xRSEit represents the value of the output of the first stage corresponding to the variable of the i-th state in the time interval t-th. The above equation is also linearized using the new inequalities (32) and (33), and the absolute nonlinear operator is removed.(32)rRSEit≥xRSEit−x^it,(33)rRSEit≥x^it−xRSEit.The 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 effects(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 first stage outputs), and consequently, improve the performance of the state estimator


## 3.4. 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 defined according to (34).(34)MSE=1NS∑i=1NSx^i−xi2.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 xi is the actual value of the i-th state variable.The 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 flow (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 flow).In this paper, according to the study of disturbances in the network at different 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 defined in the following equation (35):(35)RMSD=1NT∑t=1NT1NS∑i=1NSx^it−xit2.The RMSD value for the magnitude and phase of the voltages is calculated separately; because the scale of these two variables is different [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.


## 4. Numerical Studies

The 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 flow through the use of the MATPOWER 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. The 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. The 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 different disturbances to the network, the DSE method has been implemented in several consecutive time intervals, and in some of these time intervals, different types of disturbances have been considered. The number of time intervals is NT = 24.In examining the RSE method considered in this research (the first stage of the proposed RDSE method), the coefficientc (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, different types of scenarios are examined and analyzed to clarify the benefits of the ability to discard invalid prediction values in the event of network disturbances. Moreover, several disturbances have been assumed at different times to evaluate the performance of the proposed state estimator. In addition, different outliers have been hypothesized to investigate the effect of the robustness of the RDSE method.


### 4.1. 3-Bus Network

An overview of the 3-bus network under consideration is given in [27]. The amount of base power is assumed as Sbase=100MVA. Moreover, the production capacity of busbar 1 is equal to 0.8p.u.+0.25p.u.. Busbar 2 is assumed to be the slack busbar and busbar 1 is the PMU-equipped busbar. The 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 different time intervals. Busbar 3 load information and network line information are presented in Tables 1 and 2.Table 1
Transmission line parameters for the 3-bus network.
Transmission line numberB (p.u.)G (p.u.)X (p.u.)R (p.u.)10.02500.200.0520.02000.150.0230.02000.150.02Table 2
Busbar 3 load information in the 3-bus network.
t123456789101112131415161718192021222324P (p.u.)1.0201.0401.0611.2731.2481.2601.2861.2861.3111.2851.2591.2471.2341.2221.1980.9580.9580.9870.9770.9971.0171.0371.0581.079Q (p.u.)0.1100.1120.1150.1380.1350.1360.1390.1390.1420.1390.1360.1350.1330.1320.1290.1030.1060.1070.1060.1080.1100.1120.1140.117Due to the direct measurement of the bus angle by the PMU, the reference bus angle is not considered. The various assumptions for this network in the relevant studies are presented in Table3.Table 3
Disturbances on the 3-bus network at different times.
No.Time of disturbanceType of disturbance1t = 420% increase in the load of busbar 32t = 9Shutting down the generator of busbar 13t = 13Turning on the generator of busbar 14t = 1620% decrease in the load of busbar 35t = 20Disconnecting line 3 from the network6t = 23Inserting line 3 into the networkWhen the generator of busbar 1 is disconnected from the network, the network load and network losses are supplied by busbar 2.


### 4.2. IEEE Standard 14-Bus Network

An overview of the IEEE standard 14-bus network can be seen in [54]. Unnecessary mention of the specifications of this network (including line information, increasing costs of power plant production, transformer tap changers, and busbar information) is avoided in this research. This information is available from various sources, including reference [55].Buses 2, 6, 9, and 11 are assumed to be equipped with PMUs [54]. Therefore, 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.Table 4
Disturbances on the IEEE standard 14-bus network at different times.
No.Time of disturbanceType of disturbance1t = 5Increasing the load of busbar 132t = 12Shutting down the generator of busbar 23t = 18Disconnecting the transmission line between busbars 6 and 11 from the network


### 4.3. IEEE Standard 57-Bus Network

An overview of the IEEE standard 57-bus network can be seen in [56]. It is assumed that buses 1, 4, 6, 9, 15, 20, 24, 28, 30, 32, 36, 38, 41, 46, 50, 53, and 57 are equipped with PMUs [56] and measurements of other busbars from those obtained through the indirect measurement process. Table 5 presents the various assumptions for this network in the relevant studies.Table 5
Disturbances on the IEEE standard 57-bus network at different times.
No.Time of disturbanceType of disturbance1t = 5Increasing the load of busbars 31, 32 and 332t = 11Shutting down the generator of busbar 33t = 17Disconnecting the transmission line between busbars 18 and 19 from the network


### 4.4. Scenarios under Consideration

To evaluate the performance of the state estimators mentioned in the previous sections, several scenarios are considered in this research. The 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 algorithmScenarios 1 to 6 reflect the existing works, and scenarios 7 and 8 reflect the proposed method of this paper. Different comparisons are conducted by presenting several tables and figures. It should be noted that scenarios 5 and 6 run at one-time interval, and the rest of the scenarios run at 24-time intervals. Therefore, the RMSD index is not used to evaluate the performance of scenarios 5 and 6.


### 4.5. Numerical Studies of Optimization-Based Dynamic State Estimation Method

In this section, the performance of the state estimator under different conditions is compared. Specifically, scenarios 1 to 4 are considered in this section. The value of RMSD obtained in different cases is according to Table6.Table 6
Comparison of RMSD index values in scenarios 1 to 4.
Scenario1234AssumptionsType of SESSEDSEDSEDSEExplanationBased on optimizationBased on optimizationBased on optimizationBased on optimizationCalculation of predictions×✓✓✓Event of disturbances××✓✓Discarding the invalid predictions××✓Type of network57-bus14-bus3-bus57-bus14-bus3-bus57-bus14-bus3-bus57-bus14-bus3-busRMSDm0.00330.00740.00820.00380.00800.00910.00310.00690.00790.00420.00880.0103RMSDp0.00400.00820.00870.00420.00910.00990.00370.00810.00830.00510.00980.0126Considering 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.Table7 shows the values of binary coefficient γit for the 3-bus test network in scenario 4 at different times. The 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, final estimated values, and actual values extracted from the power flow, with an accuracy of 4 decimal places in scenario 4.Table 7
Values of binary coefficientγit for the 3-bus test network in scenario 4.
it123456789101112131415161718192021222324V1110100011111110110011011δ1110110011111110111011011V2000000000000000000000000δ2110000111000000000000001V3000000001100110110011011δ3000000001100110110011011Table 8
Comparison of measured, predicted, estimated, and true values for the voltage phasor of busbar 3 in the 3-bus network.
tTrue valuesEstimated valuesPrediction valuesMeasurement valuesV3δ3V3δ3V3δ3V3δ310.9791−0.06060.9788−0.05680.9691−0.05000.9691−0.056620.9788−0.06220.9793−0.05830.9657−0.05290.9656−0.069730.9782−0.06500.9861−0.06080.9621−0.05590.9686−0.068040.9727−0.08800.9775−0.08370.9649−0.05890.9731−0.091550.9735−0.08560.9765−0.08160.9694−0.06190.9710−0.089860.9729−0.08740.9802−0.08290.9675−0.06480.9746−0.083870.9724−0.08980.9733−0.08610.9710−0.06780.9739−0.088380.9725−0.08930.9717−0.08540.9704−0.07080.9716−0.095690.9406−0.14380.9474−0.13980.9681−0.07380.9479−0.1272100.9415−0.14150.9435−0.13720.9471−0.13050.9485−0.1250110.9423−0.13770.9459−0.13410.9450−0.13560.9494−0.1213120.9427−0.13620.9424−0.13170.9458−0.14060.9397−0.1247130.9737−0.08420.9748−0.08000.9363−0.14570.9733−0.0843140.9741−0.08280.9766−0.07910.9668−0.08870.9735−0.0782150.9746−0.08000.9810−0.07620.9699−0.09160.9732−0.0762160.9809−0.05300.9849−0.04890.9696−0.09460.9635−0.0605170.9802−0.05550.9801−0.05120.9601−0.09760.9666−0.0541180.9800−0.05630.9868−0.05230.9630−0.10060.9673−0.0514190.9803−0.05520.9843−0.05120.9637−0.10360.9663−0.0567200.9508−0.20410.9562−0.19900.9627−0.10650.9081−0.1322210.9492−0.21320.9561−0.20990.9038−0.18620.9114−0.1374220.9481−0.22000.9486−0.21630.9078−0.19120.9122−0.1401230.9783−0.06430.9794−0.06000.9086−0.19630.9675−0.0704240.9777−0.06680.9842−0.06250.9647−0.11840.9681−0.0687The following points are noteworthy in connection with Table7:(i)
The value of the binary coefficientγit corresponding to the state variable V2 (reference bus voltage magnitude) is equal to 0 at all times.(ii)
The value of the binary coefficientγit corresponding to the state variable V3 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 affected 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 coefficientγit corresponding to the state variable V3 at other time intervals is equal to 1.(iv)
The degree of influence 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.(v)
Due to the installation of PMU in busbar 1, the prediction values corresponding to the busbar 1 state variables are discarded most of the time.Tables9 and 10 show the binary coefficient values equal to 1 for PMU-free buses in the IEEE standard 14-bus and 57-bus networks in scenario 4 at different times. Similar to 3-bus networks, PMU-equipped busbars contain more accurate measurements than PMU-free busbars, and the corresponding binary coefficient is often equal to 1. Therefore, Tables 9 and 10 focus on buses without PMU.Table 9
Values of binary coefficientγit for the 14-bus test network in scenario 4.
tBusbars withγδb,t=1Busbars withγVb,t=15-69-10-11-12-139-10-11-12-1312-13-14-151-3-5—18-19-2011-12-1311-12-13Table 10
Values of binary coefficientγit for the 57-bus test network in scenario 4.
tBusbars withγδb,t=1Busbars withγVb,t=12-3-431-33-3425-31-33-34-3511-12-132-3-5-7-8-14-18-19-29-45—17-1818-1918-19Concerning Table9, the following can be mentioned:(i)
At timet = 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 affected. As expected, the busbars near the site of the disturbance were affected.(ii)
At timet = 12, bus generator 2 is tripped out. This does not constitute a sudden change in the voltage magnitude state variables, and therefore the binary coefficient γit corresponding to none is not equal to 1 at t = 12. This 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)
The outage of the line between busbars 6 and 11 at timet = 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 affects 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 specific time.Similar cases can be mentioned in connection with Table10. It is only necessary to pay attention to the fact that the occurrence of disturbances in the 57-bus network has less effect 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.


### 4.6. 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 different types of outliers. The diagrams in Figure1 show the output of the state estimation process for scenarios 5 and 6 without applying outliers.Figure 1
Comparison of true and estimated values for the magnitude and phase of the voltage of busbars in the 14-bus network by performing scenarios 5 and 6 without applying an outlier. (a) Phase of busbars voltage—scenario 5. (b) Magnitude of busbars voltage—scenario 5. (c) Phase of busbars voltage—scenario 6. (d) Magnitude of busbars voltage—scenario 6.
(a)(b)(c)(d)The diagrams in Figure2 show the output of the state estimation process for scenarios 5 and 6 by applying an outlier. The type of this piece of data is the change of the active power flow between busbars 2 and 3. When applying this outlier, it is assumed that the noise of the measurements follows the Gaussian distribution.Figure 2
Comparison of true and estimated values for the magnitude and phase of the voltage of busbars in the 14-bus network by performing scenarios 5 and 6 by applying an outlier of the type of active power change flow between busbars 2 and 3. (a) Phase of busbars voltage—scenario 5. (b) Magnitude of busbars voltage—scenario 5. (c) Phase of busbars voltage—scenario 6. (d) Magnitude of busbars voltage—scenario 6.
(a)(b)(c)(d)Figure 1 and 2 diagrams illustrate the adverse effect of the presence of outliers on the final 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 effect.Now it is assumed that the noise of the measurements is non-Gaussian, and the corresponding probability density function is obtained fromfx=1/σe−x−μ2/σ. The 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.Figure 3
Comparison of true and estimated values for the magnitude of bus voltage in the 14-bus network by performing scenarios 5 and 6 assuming Gaussian and non-Gaussian noise of measurements. (a) Magnitude of busbars voltage—scenario 5 non-Gaussian noise ME = 0.0242. (b) Magnitude of busbars voltage—scenario 5 Gaussian noise ME = 0.0111. (c) Magnitude of busbars voltage—scenario 6 non-Gaussian noise ME = 0.0109. (d) Magnitude of busbars voltage—scenario 6 Gaussian noise ME = 0.0109.
(a)(b)(c)(d)The diagrams in Figure3 show the increase in the error of the final 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 effect of non-Gaussian noise on the final output of the state estimation process is manageable.


### 4.7. 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. Table11 presents the simulation results performed on the IEEE standard 14-bus network.Table 11
Comparison of RMSD index values in scenarios 7 and 8.
Scenario78AssumptionsType of SEDSERDSEExplanationBased on the iteration algorithm and optimizationBased on the iteration algorithm and optimizationCalculation of predictions✓✓Event of disturbances✓✓Discarding the invalid predictions✓✓Type of outlier—Changing the reactive power of busbar 11 injected into the networkChanging the active power of busbar 11 injected into the networkChanging the active power flow between busbars 2 and 3—Changing the reactive power of busbar 11 injected into the networkChanging the active power of busbar 11 injected into the networkChanging the active power flow between busbars 2 and 3RMSDm0.00690.01130.00970.01010.00670.00840.00710.0074RMSDp0.00800.01210.01040.01120.00760.00910.00780.0081


#### 4.7.1. Sensitivity Analysis

The results show that the proposed RDSE method, while providing a good performance (low RMSD), can simultaneously manage the effect of network disturbances on predictions and the effect of outliers on the measurements. By examining Tables6 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. This improvement from the perspective of the RMSDm and RMSDp 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. This weakening from the perspective of the RMSDm and by assuming different 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 RMSDp and by assuming the aforementioned different 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 effect of outliers on the final output of the proposed two-stage robust dynamic state estimation process is determined. From the perspective of the RMSDm index and considering simultaneously different disturbances, a non-Gaussian noise, and aforementioned different outliers, this improvement is equal to 25.66%, 26.8%, and 26.73%, respectively. Similarly, this improvement from the perspective of the RMSDp index and considering simultaneously the aforementioned conditions is equal to 24.79%, 25%, and 27.68, respectively.Figure4 shows the output of the state estimation process for scenarios 7 and 8, assuming Gaussian and non-Gaussian noise measurements.Figure 4
Comparison of true and estimated values for the magnitude of bus voltage in the 14-bus network by performing scenarios 7 and 8 assuming Gaussian and non-Gaussian noise measurements. (a) Magnitude of busbars voltage—scenario 7 non-Gaussian noise. (b) Magnitude of busbars voltage—scenario 7 Gaussian noise. (c) Magnitude of busbars voltage—scenario 8 non-Gaussian noise. (d) Magnitude of busbars voltage—scenario 8 Gaussian noise.
(a)(b)(c)(d)Examination of recent graphs confirms the beneficial effect of using the robustness of the state estimation feature in the proposed two-stage process; in this way, the adverse effect of non-Gaussian noise on the final outputs of the two-stage process of state estimation is controlled, and the estimated values are closer to the true values.


## 4.1. 3-Bus Network

An overview of the 3-bus network under consideration is given in [27]. The amount of base power is assumed as Sbase=100MVA. Moreover, the production capacity of busbar 1 is equal to 0.8p.u.+0.25p.u.. Busbar 2 is assumed to be the slack busbar and busbar 1 is the PMU-equipped busbar. The 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 different time intervals. Busbar 3 load information and network line information are presented in Tables 1 and 2.Table 1
Transmission line parameters for the 3-bus network.
Transmission line numberB (p.u.)G (p.u.)X (p.u.)R (p.u.)10.02500.200.0520.02000.150.0230.02000.150.02Table 2
Busbar 3 load information in the 3-bus network.
t123456789101112131415161718192021222324P (p.u.)1.0201.0401.0611.2731.2481.2601.2861.2861.3111.2851.2591.2471.2341.2221.1980.9580.9580.9870.9770.9971.0171.0371.0581.079Q (p.u.)0.1100.1120.1150.1380.1350.1360.1390.1390.1420.1390.1360.1350.1330.1320.1290.1030.1060.1070.1060.1080.1100.1120.1140.117Due to the direct measurement of the bus angle by the PMU, the reference bus angle is not considered. The various assumptions for this network in the relevant studies are presented in Table3.Table 3
Disturbances on the 3-bus network at different times.
No.Time of disturbanceType of disturbance1t = 420% increase in the load of busbar 32t = 9Shutting down the generator of busbar 13t = 13Turning on the generator of busbar 14t = 1620% decrease in the load of busbar 35t = 20Disconnecting line 3 from the network6t = 23Inserting line 3 into the networkWhen the generator of busbar 1 is disconnected from the network, the network load and network losses are supplied by busbar 2.


## 4.2. IEEE Standard 14-Bus Network

An overview of the IEEE standard 14-bus network can be seen in [54]. Unnecessary mention of the specifications of this network (including line information, increasing costs of power plant production, transformer tap changers, and busbar information) is avoided in this research. This information is available from various sources, including reference [55].Buses 2, 6, 9, and 11 are assumed to be equipped with PMUs [54]. Therefore, 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.Table 4
Disturbances on the IEEE standard 14-bus network at different times.
No.Time of disturbanceType of disturbance1t = 5Increasing the load of busbar 132t = 12Shutting down the generator of busbar 23t = 18Disconnecting the transmission line between busbars 6 and 11 from the network


## 4.3. IEEE Standard 57-Bus Network

An overview of the IEEE standard 57-bus network can be seen in [56]. It is assumed that buses 1, 4, 6, 9, 15, 20, 24, 28, 30, 32, 36, 38, 41, 46, 50, 53, and 57 are equipped with PMUs [56] and measurements of other busbars from those obtained through the indirect measurement process. Table 5 presents the various assumptions for this network in the relevant studies.Table 5
Disturbances on the IEEE standard 57-bus network at different times.
No.Time of disturbanceType of disturbance1t = 5Increasing the load of busbars 31, 32 and 332t = 11Shutting down the generator of busbar 33t = 17Disconnecting the transmission line between busbars 18 and 19 from the network


## 4.4. Scenarios under Consideration

To evaluate the performance of the state estimators mentioned in the previous sections, several scenarios are considered in this research. The 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 algorithmScenarios 1 to 6 reflect the existing works, and scenarios 7 and 8 reflect the proposed method of this paper. Different comparisons are conducted by presenting several tables and figures. It should be noted that scenarios 5 and 6 run at one-time interval, and the rest of the scenarios run at 24-time intervals. Therefore, the RMSD index is not used to evaluate the performance of scenarios 5 and 6.


## 4.5. Numerical Studies of Optimization-Based Dynamic State Estimation Method

In this section, the performance of the state estimator under different conditions is compared. Specifically, scenarios 1 to 4 are considered in this section. The value of RMSD obtained in different cases is according to Table6.Table 6
Comparison of RMSD index values in scenarios 1 to 4.
Scenario1234AssumptionsType of SESSEDSEDSEDSEExplanationBased on optimizationBased on optimizationBased on optimizationBased on optimizationCalculation of predictions×✓✓✓Event of disturbances××✓✓Discarding the invalid predictions××✓Type of network57-bus14-bus3-bus57-bus14-bus3-bus57-bus14-bus3-bus57-bus14-bus3-busRMSDm0.00330.00740.00820.00380.00800.00910.00310.00690.00790.00420.00880.0103RMSDp0.00400.00820.00870.00420.00910.00990.00370.00810.00830.00510.00980.0126Considering 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.Table7 shows the values of binary coefficient γit for the 3-bus test network in scenario 4 at different times. The 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, final estimated values, and actual values extracted from the power flow, with an accuracy of 4 decimal places in scenario 4.Table 7
Values of binary coefficientγit for the 3-bus test network in scenario 4.
it123456789101112131415161718192021222324V1110100011111110110011011δ1110110011111110111011011V2000000000000000000000000δ2110000111000000000000001V3000000001100110110011011δ3000000001100110110011011Table 8
Comparison of measured, predicted, estimated, and true values for the voltage phasor of busbar 3 in the 3-bus network.
tTrue valuesEstimated valuesPrediction valuesMeasurement valuesV3δ3V3δ3V3δ3V3δ310.9791−0.06060.9788−0.05680.9691−0.05000.9691−0.056620.9788−0.06220.9793−0.05830.9657−0.05290.9656−0.069730.9782−0.06500.9861−0.06080.9621−0.05590.9686−0.068040.9727−0.08800.9775−0.08370.9649−0.05890.9731−0.091550.9735−0.08560.9765−0.08160.9694−0.06190.9710−0.089860.9729−0.08740.9802−0.08290.9675−0.06480.9746−0.083870.9724−0.08980.9733−0.08610.9710−0.06780.9739−0.088380.9725−0.08930.9717−0.08540.9704−0.07080.9716−0.095690.9406−0.14380.9474−0.13980.9681−0.07380.9479−0.1272100.9415−0.14150.9435−0.13720.9471−0.13050.9485−0.1250110.9423−0.13770.9459−0.13410.9450−0.13560.9494−0.1213120.9427−0.13620.9424−0.13170.9458−0.14060.9397−0.1247130.9737−0.08420.9748−0.08000.9363−0.14570.9733−0.0843140.9741−0.08280.9766−0.07910.9668−0.08870.9735−0.0782150.9746−0.08000.9810−0.07620.9699−0.09160.9732−0.0762160.9809−0.05300.9849−0.04890.9696−0.09460.9635−0.0605170.9802−0.05550.9801−0.05120.9601−0.09760.9666−0.0541180.9800−0.05630.9868−0.05230.9630−0.10060.9673−0.0514190.9803−0.05520.9843−0.05120.9637−0.10360.9663−0.0567200.9508−0.20410.9562−0.19900.9627−0.10650.9081−0.1322210.9492−0.21320.9561−0.20990.9038−0.18620.9114−0.1374220.9481−0.22000.9486−0.21630.9078−0.19120.9122−0.1401230.9783−0.06430.9794−0.06000.9086−0.19630.9675−0.0704240.9777−0.06680.9842−0.06250.9647−0.11840.9681−0.0687The following points are noteworthy in connection with Table7:(i)
The value of the binary coefficientγit corresponding to the state variable V2 (reference bus voltage magnitude) is equal to 0 at all times.(ii)
The value of the binary coefficientγit corresponding to the state variable V3 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 affected 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 coefficientγit corresponding to the state variable V3 at other time intervals is equal to 1.(iv)
The degree of influence 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.(v)
Due to the installation of PMU in busbar 1, the prediction values corresponding to the busbar 1 state variables are discarded most of the time.Tables9 and 10 show the binary coefficient values equal to 1 for PMU-free buses in the IEEE standard 14-bus and 57-bus networks in scenario 4 at different times. Similar to 3-bus networks, PMU-equipped busbars contain more accurate measurements than PMU-free busbars, and the corresponding binary coefficient is often equal to 1. Therefore, Tables 9 and 10 focus on buses without PMU.Table 9
Values of binary coefficientγit for the 14-bus test network in scenario 4.
tBusbars withγδb,t=1Busbars withγVb,t=15-69-10-11-12-139-10-11-12-1312-13-14-151-3-5—18-19-2011-12-1311-12-13Table 10
Values of binary coefficientγit for the 57-bus test network in scenario 4.
tBusbars withγδb,t=1Busbars withγVb,t=12-3-431-33-3425-31-33-34-3511-12-132-3-5-7-8-14-18-19-29-45—17-1818-1918-19Concerning Table9, the following can be mentioned:(i)
At timet = 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 affected. As expected, the busbars near the site of the disturbance were affected.(ii)
At timet = 12, bus generator 2 is tripped out. This does not constitute a sudden change in the voltage magnitude state variables, and therefore the binary coefficient γit corresponding to none is not equal to 1 at t = 12. This 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)
The outage of the line between busbars 6 and 11 at timet = 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 affects 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 specific time.Similar cases can be mentioned in connection with Table10. It is only necessary to pay attention to the fact that the occurrence of disturbances in the 57-bus network has less effect 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.


## 4.6. 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 different types of outliers. The diagrams in Figure1 show the output of the state estimation process for scenarios 5 and 6 without applying outliers.Figure 1
Comparison of true and estimated values for the magnitude and phase of the voltage of busbars in the 14-bus network by performing scenarios 5 and 6 without applying an outlier. (a) Phase of busbars voltage—scenario 5. (b) Magnitude of busbars voltage—scenario 5. (c) Phase of busbars voltage—scenario 6. (d) Magnitude of busbars voltage—scenario 6.
(a)(b)(c)(d)The diagrams in Figure2 show the output of the state estimation process for scenarios 5 and 6 by applying an outlier. The type of this piece of data is the change of the active power flow between busbars 2 and 3. When applying this outlier, it is assumed that the noise of the measurements follows the Gaussian distribution.Figure 2
Comparison of true and estimated values for the magnitude and phase of the voltage of busbars in the 14-bus network by performing scenarios 5 and 6 by applying an outlier of the type of active power change flow between busbars 2 and 3. (a) Phase of busbars voltage—scenario 5. (b) Magnitude of busbars voltage—scenario 5. (c) Phase of busbars voltage—scenario 6. (d) Magnitude of busbars voltage—scenario 6.
(a)(b)(c)(d)Figure 1 and 2 diagrams illustrate the adverse effect of the presence of outliers on the final 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 effect.Now it is assumed that the noise of the measurements is non-Gaussian, and the corresponding probability density function is obtained fromfx=1/σe−x−μ2/σ. The 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.Figure 3
Comparison of true and estimated values for the magnitude of bus voltage in the 14-bus network by performing scenarios 5 and 6 assuming Gaussian and non-Gaussian noise of measurements. (a) Magnitude of busbars voltage—scenario 5 non-Gaussian noise ME = 0.0242. (b) Magnitude of busbars voltage—scenario 5 Gaussian noise ME = 0.0111. (c) Magnitude of busbars voltage—scenario 6 non-Gaussian noise ME = 0.0109. (d) Magnitude of busbars voltage—scenario 6 Gaussian noise ME = 0.0109.
(a)(b)(c)(d)The diagrams in Figure3 show the increase in the error of the final 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 effect of non-Gaussian noise on the final output of the state estimation process is manageable.


## 4.7. 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. Table11 presents the simulation results performed on the IEEE standard 14-bus network.Table 11
Comparison of RMSD index values in scenarios 7 and 8.
Scenario78AssumptionsType of SEDSERDSEExplanationBased on the iteration algorithm and optimizationBased on the iteration algorithm and optimizationCalculation of predictions✓✓Event of disturbances✓✓Discarding the invalid predictions✓✓Type of outlier—Changing the reactive power of busbar 11 injected into the networkChanging the active power of busbar 11 injected into the networkChanging the active power flow between busbars 2 and 3—Changing the reactive power of busbar 11 injected into the networkChanging the active power of busbar 11 injected into the networkChanging the active power flow between busbars 2 and 3RMSDm0.00690.01130.00970.01010.00670.00840.00710.0074RMSDp0.00800.01210.01040.01120.00760.00910.00780.0081


### 4.7.1. Sensitivity Analysis

The results show that the proposed RDSE method, while providing a good performance (low RMSD), can simultaneously manage the effect of network disturbances on predictions and the effect of outliers on the measurements. By examining Tables6 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. This improvement from the perspective of the RMSDm and RMSDp 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. This weakening from the perspective of the RMSDm and by assuming different 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 RMSDp and by assuming the aforementioned different 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 effect of outliers on the final output of the proposed two-stage robust dynamic state estimation process is determined. From the perspective of the RMSDm index and considering simultaneously different disturbances, a non-Gaussian noise, and aforementioned different outliers, this improvement is equal to 25.66%, 26.8%, and 26.73%, respectively. Similarly, this improvement from the perspective of the RMSDp index and considering simultaneously the aforementioned conditions is equal to 24.79%, 25%, and 27.68, respectively.Figure4 shows the output of the state estimation process for scenarios 7 and 8, assuming Gaussian and non-Gaussian noise measurements.Figure 4
Comparison of true and estimated values for the magnitude of bus voltage in the 14-bus network by performing scenarios 7 and 8 assuming Gaussian and non-Gaussian noise measurements. (a) Magnitude of busbars voltage—scenario 7 non-Gaussian noise. (b) Magnitude of busbars voltage—scenario 7 Gaussian noise. (c) Magnitude of busbars voltage—scenario 8 non-Gaussian noise. (d) Magnitude of busbars voltage—scenario 8 Gaussian noise.
(a)(b)(c)(d)Examination of recent graphs confirms the beneficial effect of using the robustness of the state estimation feature in the proposed two-stage process; in this way, the adverse effect of non-Gaussian noise on the final outputs of the two-stage process of state estimation is controlled, and the estimated values are closer to the true values.


## 4.7.1. Sensitivity Analysis

The results show that the proposed RDSE method, while providing a good performance (low RMSD), can simultaneously manage the effect of network disturbances on predictions and the effect of outliers on the measurements. By examining Tables6 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. This improvement from the perspective of the RMSDm and RMSDp 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. This weakening from the perspective of the RMSDm and by assuming different 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 RMSDp and by assuming the aforementioned different 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 effect of outliers on the final output of the proposed two-stage robust dynamic state estimation process is determined. From the perspective of the RMSDm index and considering simultaneously different disturbances, a non-Gaussian noise, and aforementioned different outliers, this improvement is equal to 25.66%, 26.8%, and 26.73%, respectively. Similarly, this improvement from the perspective of the RMSDp index and considering simultaneously the aforementioned conditions is equal to 24.79%, 25%, and 27.68, respectively.Figure4 shows the output of the state estimation process for scenarios 7 and 8, assuming Gaussian and non-Gaussian noise measurements.Figure 4
Comparison of true and estimated values for the magnitude of bus voltage in the 14-bus network by performing scenarios 7 and 8 assuming Gaussian and non-Gaussian noise measurements. (a) Magnitude of busbars voltage—scenario 7 non-Gaussian noise. (b) Magnitude of busbars voltage—scenario 7 Gaussian noise. (c) Magnitude of busbars voltage—scenario 8 non-Gaussian noise. (d) Magnitude of busbars voltage—scenario 8 Gaussian noise.
(a)(b)(c)(d)Examination of recent graphs confirms the beneficial effect of using the robustness of the state estimation feature in the proposed two-stage process; in this way, the adverse effect of non-Gaussian noise on the final outputs of the two-stage process of state estimation is controlled, and the estimated values are closer to the true values.


## 5. 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. This 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. The effectiveness of this proposed method was verified by assuming different kinds of probable disturbances such as load fluctuations and generator outages and different types of noise (Gaussian and non-Gaussian), along with the undesirable injection of outliers into three sample networks. This method increases the accuracy and enhances the redundancy of the state estimation process. By defining the root mean square deviation (RMSD) index, which in essence and based on its definition, encompasses the definition of two other usual indices; standard deviation and mean square error, different scenarios were introduced to compare the final 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 RMSDm index was improved to the extent of 4.55%, 19.32%, and 15.91%. These three figures are for consideration of simultaneous different disturbances, a non-Gaussian noise, and different 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 RMSDp index was improved to the extent of 7.14%, 20.41%, and 17.35%, considering simultaneously the aforementioned conditions.



---

*Source: 1820478-2023-03-23.xml*