**Preface to "Optimization Methods Applied to Power Systems II"**

Electrical power systems are large and complex systems that cover the generation, transmission, distribution, and use of electrical energy. The inherent complexity of these systems is augmenting due to the increase in the generation capacity from renewable energy sources. Therefore, the analysis, design, and operation of current and future electrical power systems require an efficient approach to different optimization problems, including load flow, fault location, contingency analysis, system restoration after blackouts, economic dispatch, unit commitment, islanding detection, filter design, and smart grids control. The fast evolution of this field and the difficulty of managing large systems demand the design and validation of advanced techniques for supporting decision-making processes. Currently, power system analysis, design, operation and management tasks are being supported by advances in computational optimization methods, from classical optimization techniques such as linear and nonlinear programming and integer and mixed-integer programming to recent advances in meta-heuristic optimization, including bio-inspired meta-heuristics which have allowed researchers and practitioners to obtain quality solutions in reduced response times. Taking into account the high dimensions and complexity of real-life power systems, efficient network planning, operation, or maintenance requires a permanent research effort.

Therefore, this book includes recent advances in the application of optimization techniques that directly apply to electrical power systems so that readers may familiarize themselves with new methodologies directly explained by experts in this scientific field.

> **Francisco G. Montoya, Ra ´ul Ba ˜nos Navarro** *Editors*

## *Article* **Online Economic Re-dispatch to Mitigate Line Overloads after Line and Generation Contingencies**

#### **Oswaldo Arenas-Crespo 1, John E. Candelo-Becerra 1,\* and Fredy E. Hoyos Velasco <sup>2</sup>**


Received: 15 February 2019; Accepted: 9 March 2019; Published: 13 March 2019

**Abstract:** This paper presents an online economic re-dispatch scheme based on the generation cost optimization with security constraints to mitigate line overloads before and after line and generation contingencies. The proposed optimization model considers simplification of mathematical expressions calculated from online variables as the power transfer distribution factor (PTDF) and line outage distribution factor (LODF). Thus, a first algorithm that calculates economic re-dispatch for online operation to avoid overloads during the normal operation and a second algorithm that calculates online emergency economic re-dispatch when overloads occur due to line and generator contingencies are proposed in this paper. The results show that the proposed algorithms avoid overload before and after contingencies, improving power system security, and at the same time reducing operational costs. This scheme allows a reduction of power generation units in the electricity market during online operation that considers line overloads in the power system.

**Keywords:** electricity markets; economic dispatch; optimization; maximum capacity; online operation; power system; unit commitment

#### **1. Introduction**

Power system security is an important subject for electricity companies that operate transmission networks, which depend on the system operating condition and contingencies and can lead to interruption of customer service [1–3]. Thus, some methodologies seek to identify the safe operation of the power system [2]. Although these procedures allow a safer operation, the large number of exogenous variables, such as maintenance that exceeds the planned time, emergency outage of transmission equipment, and volatility of renewable sources, lead the operation to a risk zone and new emergency actions.

A classic economic dispatch is commonly used to identify the generation costs of each unit, which considers constraints as Kirchhoff's first and second laws with AC or DC models [4]. The security-constrained economic dispatch (SCED) could include physical limits of equipment, operative limits, and contingencies [5], and the solution is found with a security-constrained optimal power flow (SCOPF), in which the objective function can be the minimization of losses, minimization of power demand rationing, minimization of the cost of operation, and others. Security constraints imply the modeling of the power grid, which leads to an important high computational cost [6]; therefore, economic dispatch models with security constraints are usually implemented by means of metaheuristic methods that imply lengthy execution time, or by means of exact models such as the interior point algorithm, which sometimes causes difficulties with convergence [7].

With this, the economic re-dispatch is made to avoid overloads at the lowest cost, as implemented in [8,9]. In this case, particle swarm optimization (PSO) is used to optimize as a stable algorithm with respect to the convergence of the non-linear modeling of the power flow equations.

Other authors have implemented a similar economic dispatch model with classic constraints by using learning machines [10]. In addition, an interesting algorithm was presented in [11], which selects the generation to participate in an economic dispatch using a direct acyclic graph (DAG). This model has been proposed as an alternative for large networks and various operational areas. In [12], a model that considers mixed integer linear programming (MILP) to minimize switch opening as a solution to reduce overloads is implemented. Further, in [13] the overloads are reduced through fuzzy logic, whose model tries to recreate the actions of the network operators; however, this model does not take into account the generation cost.

In [14], the online economic dispatch is implemented using the metaheuristics of Fuzzy Logic and Genetic Algorithm, which avoids the modeling of the complete AC system and the problems of non-convergence of exact solution methods. On the other hand, an exact solution method based on Primal Dual interior Point (PDIP) can be used [15]; however, in order to execute the final solution algorithm, different simplification stages of the dispatch scenarios to be executed by clustering can be used, so that the appropriate selection of the scenarios will have a direct impact on the consolidation of the global optimum and compliance with the system constraints.

Middle- and long-term solutions seek to guarantee security and reliability within the network. On the one hand, middle-term solutions use contingency analysis through simulations to find network expansion that mitigate the constraints found in the analysis. On the other hand, short-term solutions and online operation are based on controlling constraints where the operator or intelligent algorithms monitor and control system variables. However, during online operation there is no time to perform network expansion; here, techniques consider variation of the network topology, new economic dispatches, stress control, load shedding, and others.

Security-constrained economic dispatch (SCED) uses algorithms that are of high computational cost, so they are not normally used in the online operation. Therefore, despite previous plans for the operation through middle- and short-term analysis to carry out the secure operation of power systems, deviations of the programmed resources and the projected variables are evidenced in the operation. Therefore, this paper proposes an online economic re-dispatch to mitigate overloads of transmission elements after N−1contingencies to reduce the risk of collapse [16]. Because of the number of variables that contain the problem, the mathematical formulation can be simplified by distribution factors (DFs) such as the power transfer distribution factor (PTDF), the line outage distribution factor (LODF) [17], and the outage transfer distribution factors (OTDF) [4,5,18–20]. This technique has an advantage that allows the linearization of the power flow's equations around an operation point. Thus, the contributions of the paper are defined as follows:


The rest of the paper is organized in three more sections: Section 2 describes the main mathematical model of the online economic re-dispatch for the normal operation and the model of the online economic re-dispatch proposed for overloads; the sensitivity factors, which are applicable to the power flow and the online operation are included in the formulation. Section 3 presents the power system test case and the results obtained by applying the formulation. Finally, Section 4 includes the conclusions, recommendations and the future work in this research.

#### **2. Methodology**

In this section the mathematical formulation used for the application of the economic re-dispatch schemes, and the simplification of the formulation to apply the electrical constraints of the system are presented. In addition, we present the algorithms used for the economic re-dispatch in normal operation and the economic re-dispatch after contingencies.

#### *2.1. Economic Re-Dispatch with Power Demand Rationing Cost*

The online economic dispatch model to reduce overloads after N−1contingencies can be formulated as in Equation (1), where Δ*PGi* is the power generation change in bus *i*, *CGi* is the operational generation cost in bus *i*, Δ*PRi* power load shedding change in bus *i*, *CRi* is the load shedding cost in bus *i*, and *i* = 1, 2, 3 ... *m* buses. This equation corresponds to the objective function, which is used to minimize the operating cost of the system. The expression considers the sum of the generation cost and the cost of power demand rationing. It should be noted that the generation cost (*CGi*) is assumed as a constant value that is expressed in terms of the cost per power supplied (\$/MWh), as well as the cost of power demand rationing. However, the latter is much greater than the generation cost and it can be considered as a constant value of great magnitude, because it is only required to help the convergence of the model:

$$Minimize\ F = \sum\_{i=1}^{m} (\Delta PG\_i \* CG\_i + \Delta PR\_i \* CR\_i). \tag{1}$$

The model starts with an initial condition, dispatching the central generation plants (*PG*<sup>0</sup> *i* ), the power flowing through transmission elements (*P*<sup>0</sup> *<sup>i</sup>* ) and the power demanded in the load buses of the power system (*PDi*). In the present work, it is assumed that all these initial parameters are obtained online by measurement units as the supervisory control and data acquisition (SCADA), because this information can be managed from this type of application.

When there is no feasible safe dispatch and the power demand rationing can be applied to bus *i*, the model would reflect that situation, obtaining a value of Δ*PRi* greater than zero. In this case, the final power demanded in bus *i* (*PD i* ) would be reduced as expressed in Equation (2):

$$PD\_i' = \ \_PD\_i - \Delta PR\_i. \tag{2}$$

Therefore, the decision variable for power rationing (Δ*PRi*) is modeled as a generator in each bus with power limits on the demand at each bus, as expressed in Equation (3):

$$0 \le \Delta PR\_i \le PD\_{i\prime} \text{ i } = 1, 2, \dots, m. \tag{3}$$

The active generation or generation synchronized to the system is the only one that should be considered in the optimization model, because it is the only one that can be managed in real time. This generation could consider renewable generation sources such as solar, wind, and even battery banks (BESS).

Starting from the initial operation (*PG*<sup>0</sup> *<sup>i</sup>* , *<sup>P</sup>*<sup>0</sup> *<sup>i</sup>* , *and PDi*), the model calculates the deltas of generation (Δ*PGi*), which are required to maintain a preventive condition where the N−1contingencies do not generate overloads. Once the model converges, the new power for each generator in bus *i* (Δ*PGi*) will have a value calculated by Equation (4):

$$
\Delta PG\_l t = \, P G\_l^0 + \Delta P G\_l. \tag{4}
$$

Equation (5) is used to prevent the total generation deviation from exceeding the allowed level of automatic generation control (AGC) for the online re-dispatch algorithm, where *ResAGC* represents the power reserves used for generation re-dispatch. This value of AGC can be managed by the *ResAGC* parameter, which could even be 0 MW. Equation (6) is used to avoid the operating limits of each generator are violated, where *PGmaxi* and *PGmini* are the maximum and minimum power generation in bus *i*, respectively:

$$-\operatorname{Res} A \mathbf{G} \mathbf{C} \le \sum\_{k=1}^{m} \Lambda P \mathbf{G}\_{i} \le \operatorname{Res} A \mathbf{G} \mathbf{C}, \ i = 1, 2, \dots, m. \tag{5}$$

$$PGmin\_{l} \leq PG\_{l}^{0} + \Delta PG\_{l} \leq PGmax\_{l}, i = 1, 2, \dots, m. \tag{6}$$

The power that flows for each transmission element (*Pi*) is calculated from the initial state *P*<sup>0</sup> *i* , which corresponds to the power measured online at transmission lines and transformers. To this value, the calculated delta of the power flow (Δ*PGk*) and delta power rationing (Δ*PRk*) are added, multiplied by the power transfer distribution factors PTDF that relate to the contribution of the power injection in the bus with respect to the transmission element (*λik*), as shown in Equation (7):

$$P\_i = P\_i^0 + \sum\_{k=1}^m \lambda\_{ik} \* \Delta P G\_k + \sum\_{k=1}^m \lambda\_{ik} \* \Delta P R\_k \; , i = 1, 2, \dots, n. \tag{7}$$

Equation (8) corresponds to the security constraints that prevent the transmission elements from being overloaded above their allowed values, after generation changes and without contingencies; where *Pnomi* is the rated power flow of each transmission element *i* without contingencies:

$$-\operatorname{Phom}\_{i} \le \operatorname{P}\_{i} \le \operatorname{Phom}\_{i}, i = 1, 2\dots, n. \tag{8}$$

Finally, the constraints in Equation (9), ensure that no N−1contingency leads to the overload of other elements, by including the sensitivity factor LODF (*γij*) between two lines. Herein, the term *Pmax* refers to the maximum power flow for transmission element *i* under single contingencies:

$$\begin{array}{l} -P\text{max}\_{i} \le |P\_{i} + P\_{j} \ast \gamma\_{ij}| \le \text{ } P\text{max}\_{i} \\\ j = 1, 2, \ldots, n, \ j = 1, 2, \ldots, n, \ i \ne j. \end{array} \tag{9}$$

#### *2.2. Economic Re-Dispatch without Power Demand Rationing Cost*

If an economic dispatch solution is required without power rationing, the objective function is defined as presented in Equation (10). This function considers minimizing the operating cost of the system:

$$Minimize\ F = \sum\_{i=1}^{m} (\Delta P G\_i \* \mathbb{C}G\_i) + \sum\_{i=1}^{n} \alpha \* PA d\_i. \tag{10}$$

The constraints defined above are maintained, and the power flowing through each transmission element is calculated from the initial power of the system, as shown in Equation (11):

$$P\_i = P\_i^0 + \sum\_{k=1}^m \lambda\_{ik} \* \Delta P G\_k \; , i = 1, 2, \dots, n. \tag{11}$$

In addition, Equation (2) is delete and Equation (3) is modified with the expression shown in Equation (12):

$$\begin{array}{c} -(Pmax\_i + PAd\_i) \le |P\_i + P\_j \* \gamma\_{ij}| \le |Pmax\_i + PAd\_i|\\ \text{i} = 1, 2, \dots, n, \text{ j} = 1, 2, \dots, n, \text{ i } \ne j. \end{array} \tag{12}$$

With this modification the model that considers electricity rationing is eliminated, and now the capacity of the transmission elements after contingencies is considered in the formulation. The variable *PAdi* allows the convergence of the model, even when the system is unsecure with critical contingencies, allowing security to be maximized. The critical contingencies correspond to all those contingencies that generate overloads and that said overload cannot be alleviated by any feasible online economic dispatch. Therefore, the result of this model corresponds to the most secure economic dispatch possible without considering the rationing and improving the capacity of the transmission elements to contingency.

Because the process is iterative according to the total number of critical contingencies, this implies that the model must be linearized around the point of operation whenever the topology of the system is altered from the evaluation of each critical contingency. The above implies that the matrix of PTDF sensitivity factors must be updated for each topology. However, there is another possibility to update the sensitivity matrix PTDF, by replacing this matrix with the calculation of the sensitivities known as OTDF [21].

According to this, the power flow for transmission elements that was considered in Equation (7) is modified and the new expression is shown in Equation (13), where:

$$P\_i = P\_i^0 + \sum\_{k=1}^m \psi\_{ik} \* \Delta P G\_k + \sum\_{k=1}^m \psi\_{ik} \* \Delta P R\_k \; , i = 1, 2, \dots, n. \tag{13}$$

where *ψik* represents the sensitivity factor OTDF, that models the power flow distribution in a bus *k*, to the monitored element *i*, considering the outage of an element *j* (for our case the critical contingency evaluated). Therefore, the calculation of the sensitivities for the modelling variations in the power flow after various topological changes becomes simpler, if it is formulated from the OTDF sensitivities matrix instead of PTDF.

#### *2.3. Algorithm for the Online Economic Re-Dispatch*

The general model of the online economic re-dispatch algorithm used in the research is presented in Figure 1. The interaction between the power system information SCADA, the AC simulator, and the linear optimization model is observed.

**Figure 1.** First algorithm for the online economic dispatch.

The model is divided into three main stages. The first one (inputs) corresponds to the capture of information corresponding to the physical system. The second one (process) corresponds to the evaluation of this information and integration of the optimization models and AC simulator. Finally, the third one (outputs) corresponds to the economic re-dispatch solution in each operating state which could be implemented in an automated system or through operational actions of operators.

In the first stage of the complete methodology, the inputs of the system are considered as the required variables for the procedure. The input E1 corresponds to the status of the transmission lines and transformers, representing the elements that are in service (ON = 1) or out of service (OFF = 0). The input E2 defines the difference between the forecasted demand and the actual demand, representing the power deviation of the system. The input E3 corresponds to the difference between the scheduled power dispatch and the current power dispatch presented in the online transaction; for this, each generator must be monitored to calculate the deviation. The input E4 represents the price of each generator for the online operation period. The input E5 refers to the minimum and maximum limits of the generation plants, which, in the online operation, are considered variables, mainly in thermal plants and will depend on factors such as generator temperature, generator's configuration (cycle center combined), start ramps, etc. The input E6 is the linearized model obtained from the current power system condition, and refers to the values of real power that flows through the various transmission elements that are monitored; it is not necessary to monitor or model the complete power system for use in the economic re-dispatch with overloads and after contingencies.

The second stage corresponds to the logic and subroutines of the general process. For example, the process L1 identifies the state of a line or transformer defined as in service (ON = 1) and out of service (OFF = 0). This logic can be implemented from measures of each element or switch positions. The process L2 corresponds to the logic in which the presence of a deviation of generation (E1) or demand (E2) is evaluated, so that the algorithm recognizes the presence of a deviation, it must be fulfilled; thus, if a deviation is greater than a value *ArP*, then the logic L1 can be expressed as the inequality Δ*Pi* > *ArrP*.

Now, P1 refers to the economic re-dispatch optimization model with constraints and N−1contingencies. Herein, the optimization model does not consider the calculation of power demand rationing. Therefore, the optimization model defined in P1 uses the objective function of Equation (10) and the constraints of Equations (11), (12), and (13). P1 also corresponds to the main subprocess of the algorithm, because the solution of the economic re-dispatch is obtained from it, which is refined through the iterations with the AC simulator. Then, the process P2 includes the sensitivity factors PTDF and LODF, and the optimization model P1 is linearized from a current operating point. This process must be carried out whenever the topology of the network is modified, due to the deviations in the distribution of the power flows.

In the process P3 the operating state of the AC simulator must be coordinated, which implies that the state of lines and transformers in the simulator must be the same as in the real power system, as well as the power demand in the bus of the power system. However, the power dispatch to be implemented in the simulator must be that obtained through the process optimization model P1. Once the simulator is coordinated, N−1contingencies are executed with AC load flow. From the N−1contingency results obtained from process P3, the process L3 verifies that no overloads are presented in the elements. If an overload is detected in this process, then the process P4 is executed, correcting the maximum power overload with this change in the parameter of the optimization model; the process P1 is executed again thus obtaining an economic re-dispatch result, starting from the new iteration. If no overloads are detected in the process L3, then the final value is obtained.

#### *2.4. Online Emergency Economic Re-Dispatch Algorithm*

The model presented in Figure 1 should not generate power demand rationing as a solution to the problem. Instead, an adaptable optimization model is necessary for the calculation of the secure dispatch and lowest cost without sacrificing the power demand and minimizing the overload of transmission elements after N−1contingencies. However, if critical contingencies are detected when applying the algorithm presented in Figure 1, then the complementary algorithm presented in Figure 2 is executed. This new processes perform the calculation of a new dispatch and emergency load shedding, for each one of the critical contingencies found in the P1 process. That is, in this second phase of the algorithm, a new dispatch is calculated with possible load shedding for each critical contingency. This complementary sub-algorithm for critical contingencies corresponds to an optional algorithm, and is useful in cases in which the network is not fully covered before all the possible N−1contingencies, or is useful in cases in which the only solution to avoid collapse is the load shedding.

According to the above, it is assumed that an emergency load shedding is considered because of critical contingencies. The previous algorithm could have an operation similar to that of a systemic protection supplementary scheme (SPS), but considering the minimization of the generation cost and minimization of load shedding, depending on the various operation points and topologies presented in online operation. The power demand rationing should always be the minimum possible. For this the complementary algorithm of Figure 2 has been proposed, which, given the initial dispatch calculated in process P1, obtains the value of the new dispatch and emergency load shedding for each critical contingency.

**Figure 2.** Second algorithm for the online dispatch with load shedding after critical contingencies.

In the process L4, the presence of critical contingencies is detected through the evaluation of the variable (*PAdi*). If there is no risk due to critical contingencies, the algorithm behaves as shown in Figure 1; otherwise, the complementary algorithm shown in Figure 2, composed of P5 and P6, is executed. With the process P5, a most secure generation dispatch is obtained, despite the existence of critical contingencies. Additionally, the process L3 and processes P3 and P4—through which the economic re-dispatch is obtained from P1—is refined by iteration with the AC simulator.

Then, the process P6 considers an adapted optimization model, this time, with the objective of providing a solution of an economic re-dispatch with the possibility of load shedding for each of the critical contingencies. Therefore, the objective function formulated by Equation (1) is used, considering the constraints of the model presented from Equations (2) to (6). Additionally, the constraint in Equation (7) is not considered in the first part of the algorithm, because the calculation of the new economic re-dispatch and load shedding is performed after detecting critical contingencies. Thus, the optimization model is applied to the degraded network after critical contingency; therefore, from each critical contingency, a result of the new economic dispatch and different load shedding is obtained. Finally, the constraint presented in Equation (10) is not required for process P6.

#### **3. Results**

#### *3.1. IEEE 39–Bus Power System Test Case*

Figure 3 shows the IEEE 39–bus power system test case, which is a simplified model of the New England power grid. This system has ten generators; generator one is an equivalent representation of the rest of the power grid. This power system topology allows one to assess multiple power demand and generation dispatch scenarios.

**Figure 3.** IEEE 39–bus power system test case.

#### *3.2. Input Parameters*

Table 1 presents the demand value and initial economic dispatch considered in the research. This dispatch has been calculated using the optimization model with the aim of starting from a secure point, which will be altered before topology changes and generator outages. This is also done to obtain the economic re-dispatch using the proposed online economic dispatch algorithm.

Table 1 shows that the power demand rationing cost tends to an infinite value, and the optimization must search to minimize this value modifying each bus of the system. For example, in Colombia the power demand rationing cost oscillates between 130 and 2608 USD/MWh [22]. From the generation dispatch and demand presented in Table 1, the initial operating state of the power system is obtained as shown in Figure 4, in which each line and transformer obtains a color according to the level of overload presented between 0% and 100% of its maximum loadability. In this case the load flow is simulated, and the power system is represented under normal operating conditions and considering all the elements in service.

By representing the maximum loadability of transmission elements, considering N−1contingencies, and using the same method of coloring elements used in Figure 4, the heat diagram of Figure 5 is obtained, in which branches increased overloading risks when N−1contingencies occurred. Additionally, when comparing Figures 4 and 5, loadability has increased with N−1contingencies and that affects power system security. For example, line 23–24, under normal operating conditions has an overload of less than 50% and the power flow can increase to approximately twice the current value. However, for the same line after N−1contingencies, loadability reaches values higher than 90%, implying a critical state.


**Table 1.** Initial operation of the power system.

**Figure 4.** Maximum loadability without N−1contingencies.

**Figure 5.** Maximum loadability with N−1contingencies.

As discussed above, the optimization model uses security constraints based on sensitivity factors, such as: PTDF, LODF and OTDF. The main advantages of using these sensitivity factors are: first, they simplify the security constraint model for both power flow and contingencies calculations; and second, they include small errors in the optimization result when the system works with AC model. In this method, the constraint calculation is made based on the initial and real states of the power system, which means that a linearized model is achieved around the operating point and a simple model is obtained to solve the difficult problem of the economic dispatch with online security constraints. For example, Figure 6 shows a comparison between the AC and DC power flow for a 230 kV line with 260 km, obtaining that the largest error of the calculation is given by the angular opening and the voltage along the line. Therefore, the error when using sensitivity factors is small, because the problem is linearized around the real state of the power system.

**Figure 6.** Comparison between the DC power flow and the AC power flow with sensitivity factors with three different voltage magnitudes.

#### *3.3. Line Outages*

The algorithm must again perform a generation dispatch calculation based on the changes presented in the power system (L1 and L2 of Figure 1). In the case of the unexpected line outage, logic L1 will activate the algorithm, which will obtain a new economic re-dispatch. In accordance with the above, permanent line outage events that do not represent critical contingencies activate the first economic dispatch algorithm presented in Figure 1 and do not activate the second economic dispatch algorithm presented in Figure 2. For example, in the IEEE 39–bus power system the lines outages are: (a) Line 3–18, (b) Line 16–21, and (c) Line 17–18, and those allow calculating an economic re-dispatch.

For each event a, b and c, the algorithm calculates the economic dispatch and the results are shown in Table 2. The economic dispatch before the outage of line 16–21 is approximately equal to the dispatch of the base case. Therefore, it is not necessary to implement a generation re-dispatch for this case. On the other hand, the output of lines 3–18 and 17–18 originate a more relevant change in the generation dispatch. Additionally, generation plants G4, G5, G6, and G7 remain almost invariant before the outage of these three lines. As generators G4, G5, and G6 are the most economic units, then those generators are not dispatched again because of security constraints related to overloads.


**Table 2.** Results of the economic re-dispatch after the outage of lines 03–18, 16–21, and 17–18.

None of the three outages of lines 3–16, 16–21, and 17–18 create a risk situation after N−1contingencies. However, the resulting economic re-dispatch corresponds to an opportunity to save generation costs, because each outage of each of these lines creates a different topology, so that power system security changes. Thus, the algorithm performs online monitoring of the cost reductions, avoiding the overload risk after N−1contingencies, and contributing to the improvement of power system security and operating costs. By considering hour-operation periods, the outage of line 17–18 saves around USD 12,316; the outage of line 16–21 saves only USD 42.15; and the outage of line 03–18 saves USD 8859.73. According to the previous results, we can say that before the unexpected outage of lines 03–18 and 17–18, it is not necessary to carry out a re-dispatch if only the security constraints after N−1contingencies are required. However, the resulting economic re-dispatch corresponds to an opportunity to save the operating cost.

Finally, the verification of the deviation value of each of the new economic re-dispatch is important because the power delta lost or gained after each dispatch could affect the operation of the secondary frequency control or AGC. In addition, the generation delta should not be very large, since the AGC must have available a power reserve to respond to frequency events, and in this case the permissive availability of AGC of + −10 MW has been assumed. For the outage of line 03–18 there is a power generation deviation of −4.857 MW, for the outage of line 16–21 there is a deviation of 0.688 MW, and for the outage of line 17–18 there is a deviation of −2,397 MW.

#### *3.4. Critical Contingencies*

To involve the operation of the algorithm presented in Figure 2, two lines have been selected, which correspond to the outages of lines 05–08 and 13–14, and which cause critical contingencies N−1. Therefore, the algorithm will present two operation stages: the first is presented in Figure 1, where the algorithm must calculate the secure operation and most economical re-dispatch without considering load shedding, and without minimizing the overload after N−1contingency. In the second stage presented in Figure 2, the algorithm must calculate the new redistribution and load shedding corresponding to each of the critical contingencies that remain unprotected in the first stage.

In accordance with the above, Figure 7 shows the results of the economic re-dispatch of the first stage, which allows the lowest possible overload after N−1contingencies at the lowest generation cost. However, as stated above, the evaluated line outages create critical contingencies N−1. Thus, before the outage of any line (05–08 and 13–14), there is no totally secure economic re-dispatch; that is, under this condition there will be at least one N−1contingency that causes an overload above the maximum limit of the overloaded element. In Table 3, we present in detail the values of the economic re-dispatch to consider, before the outage of lines 05–08 or 13–14.


**Table 3.** Results of the economic re-dispatch after outage of lines 05-08 and 13–14.

In relation to the deviations of the generation cost, when the outage of line 13–14 is presented, the cost increases by USD 3271.93; therefore, in this case, the algorithm detects the need to increase the generation of more expensive security to reduce the risk of contingencies N−1. On the other hand, the result of the cost of generation before the exit of the line 05–08 leads to a reduction of the same in USD 9818.76, which represents an operating benefit, and at the same time leads to a reduction of risk due to overload between N−1contingencies.

When evaluating N−1contingencies before the outage of line 05–08 considering the economic re-dispatch presented in Table 3, it is observed that lines 09–39, 08–09 and 06–07 present an overload on their maximum values, as seen in Figure 7.

After the outage of line 05–08, the algorithm calculates a fictitious increase in the capacity of lines 09–39, 08–09 and 06–07, resulting in 42.3 MW, 28.2 MW and 142.5 MW, respectively. The fictitious increase in the capacity of transmission lines is calculated by means of the variable *PAdi*, of Equations (8) and (10). From the process P5 of Figure 2, critical contingencies are obtained when considering the outage of line 05–08. For this case, the critical contingencies correspond to the lines: 08–09, 06–07, 09–39 and 01–39, as shown in Figure 8.

**Figure 7.** Simulation of maximum loadability of elements when considering N−1contingencies after the outage of line 05-08.

**Figure 8.** Emergency dispatch after a contingency that considers the outage of line 05–08.

Table 4 shows the results obtained by simulating critical N−1contingencies, considering the outage of line 05–08. The loadability of each element of the power system is presented. For each critical contingency, it is required to carry out the dispatch presented in Figure 9, as well as the corresponding load shedding as shown in Table 4.


As shown in Table 4, the contingency of line 01–39, would be the only critical contingency that would not increase the operation cost. The other contingencies cause a considerable increase in the operation cost, due to the fact that they require load shedding.

#### *3.5. Economic Re-Dispatch Algorithm Operating after Generation Tripping*

For this case, a generation tripping event is considered, which activates the online economic re-dispatch algorithm using the L2 logic, presented in Figure 1, for this case a complete outage of generator G8 in bus 37. Figure 9 shows the loadability status of the transmission elements after the outage of G8 and with N−1contingencies. High loadability are presented in lines 01–39, 01–02, 02–25, 09–39, 08–09, because the slack bus, which represents the generator with AGC control of the power system, responds to the loss of the G8 generation with a real power increase and generates power congestion on those transmission lines. With an event that creates generation-demand imbalance, as the outage of G8, the redistribution of the AGC generators is required, greater than the value considered in the *ResAGC* parameter of Equation (5). Because after the outage of a non-radial transmission line, there would be no generation-demand imbalance; therefore, the adjustment value of the AGC (*ResAGC*) may be lower than after the event that generates a generation-demand imbalance.

For the event of the outage of G8, the *ResAGC* value has been modified by 248.729 MW, which corresponds to the value of the G8 generator dispatch after the outage of G8. Figure 10 shows the variation of the economic re-dispatch calculated from the generation tripping (G8), with which the system will be safe again with N−1contingency. Table 5 shows the values of the economic re-dispatch calculated after the outage of G8, a generation deviation of −19.58 MW is presented, which means that an increase of the *ResAGC* parameter is required, given that the most economical solution is the generators G2 and G3, generating an unbalance greater than 10 MW.

It should be noted that for the simplicity of this model, a cost of 0.0 USD/MWh has been assumed for the AGC generator, because the generator G1 represents not only the AGC, but also the rest of the equivalent system; therefore, the AGC is assumed to be part of the equivalent system and for simplicity the cost of this equivalent generator is assumed to be zero. Deviation of the total generation is −19.58 MW and the generation cost is USD −4945.

**Figure 9.** Maximum loadability of elements after outage of G8.

**Figure 10.** Maximum loadability of elements after outage of G8, with the new economic re-dispatch algorithm.

\*a

**Table 5.** Economic re-dispatch values after the outage of generator G8.


#### **4. Conclusions**

This paper presented an online generation dispatch scheme to minimize generation costs and the loadability of elements in the power system. The proposed method incorporated the power sensitivity factors (LODF and PTDF) in an optimization model with variables obtained from online operation, considering the generation cost, and with overload constraints. The results showed that the online economic re-dispatch identifies the operating cost reduction and reduces the overload risk created by N−1contingencies, contributing to a safety operation, and at the same time reducing the operating costs of the power system. In addition, it is a useful proposal when there are highly variable power plants such as power systems with solar and wind power.

This operating scheme could also be very useful in power systems in which there are usually appreciable differences in demand, dispatch or topology, between the scheduled operation and the actual operation. The complementary algorithm to deal with critical N−1contingencies, corresponds to a useful proposal in cases in which the network is not fully covered before all the possible N−1contingencies. It is also considered a necessary operation scheme, in cases in which the only solution to avoid collapse is load shedding. Therefore, instead of executing a programmed power demand rationing, the model uses a prediction of the required load shedding and emergency economic re-dispatch, which would be implemented only in cases of critical contingency.

The model considered sensitivity factors in order to simplify calculations, and it provided an advantage when considering large-scale systems divided into many operational areas. The online economic re-dispatch presented in the present work does not require metaheuristic techniques, so the minimization of the generation cost is done by an exact and simple method, which means it is suited to an online operation. Thanks to the model based on the objective function with the variable *PADi*, it is possible to achieve the convergence of the model, even when the system is unsafe with critical contingencies. Therefore, the result of this model corresponds to the security constraint generation economic re-dispatch without the need to calculate power demand rationing.

**Author Contributions:** Conceptualization, methodology, software, validation, formal analysis, investigation, and writing original draft were performed by O.A.-C. and J.E.C.-B. Formal analysis, writing, review, and editing was performed by F.E.H.V.

**Funding:** This research received no external funding.

**Acknowledgments:** This work was supported by the Universidad Nacional de Colombia, Sede Medellín. We would like to thank to the Department of Electrical Engineering and Automation and the Grupo de Investigación en Tecnologías Aplicadas GITA for the continuous support to our research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **An Expeditious Methodology to Assess the Effects of Intermittent Generation on Power Systems**

#### **Gracita Batista Rosas 1, Elizete Maria Lourenço 2,\*, Djalma Mosqueira Falcão <sup>3</sup> and Thelma Solange Piazza Fernandes <sup>2</sup>**


Received: 2 March 2019; Accepted: 21 March 2019; Published: 23 March 2019

**Abstract:** This paper proposes an expeditious methodology that provides hourly assessments of the effect of intermittent wind and solar power generation on the electrical quantities characterizing power systems. Currents are measured via circuit breakers to confirm the correct sizing of devices based on their rated currents. Nodal voltage magnitudes are assessed for compliance with limits imposed by regulatory authorities, whereas the active power produced by hydroelectrical generators is assessed for reserve energy. The proposed methodology leverages a fuzzy extended deterministic optimal power flow that uses in power balance equations the average hourly values of active power generated by wind and solar sources as well as hourly energy load. The power grid is modeled at the substation level to directly obtain power flow through circuit breakers. Uncertainties in power system electrical quantities are assessed for an optimal solution using a Taylor series associated with deviations from the average values of the active power produced by the wind and solar sources. These deviations are represented using a fuzzy triangular model reflecting the approximations of the probability density functions of these powers. The methodology takes into account a subjective investigation that focuses on the qualitative characteristic of these energy sources' behaviors.

**Keywords:** circuit breaker; network modeling; optimal power flow; wind and solar energies

#### **1. Introduction**

The challenges of the expansion and operation of Electrical Power Systems (EPS) have intensified over time owing to the incessant increase in the complexity of these systems. Nowadays, additional factors have been incorporated in these challenges, including uncertainties in the generation and consumption of electric energy as dictated by the growing participation of intermittent sources of energy, domestic electric vehicles, and energy storage systems. Because of these uncertainties, the technical specifications of the EPS components must be assessed, especially those associated with substations' internal devices, as is the case with circuit breakers. Therefore, it is of fundamental importance that the tools in development that seek to model these uncertainties are also capable of assessing the loading of these devices in their analysis processes.

Optimal power flow (OPF) is recognized as a powerful tool for planning an operation and has gained importance in analysis involving the incorporation of distributed generation [1–3]. OPF is able to indicate the ideal size of these generation units to avoid compromising voltage stability [4,5]. The uncertainties of wind and solar energy sources have been widely explored in OPFs, whether these were coordinated with conventional generation or in association with energy storage systems [6–8]. In this context, some works have explored the characteristics of the OPF to analyze the impact of the expansion of distributed generation on transmission line loading [9], on the presence of domestic electric vehicles in power grids [10,11], on high-voltage direct current connections [12,13], and in frequency control actions [14].

The uncertainties of energy sources are often represented in an OPF by probabilistic approaches, which are unable to incorporate qualitative information residing in human knowledge. In this case, fuzzy logic has been used to represent qualitative, vague, or incomplete knowledge in mathematical models in several areas of knowledge, including the OPF [15]. The literature presents the use of fuzzy logic combined with OPF as promising for EPS analysis [16], mainly in the definition of functions with conflicting objectives [17,18], in operational limit constraints [19,20], and in solutions that define the date of need, location, and sizing of new energy resources [21]. Some techniques incorporate fuzzy logic in the OPF solution through a linear relation [22], or through sensitivity analyses [23]. This allows for the assessment of the influence of uncertainties in energy generation in the electrical quantities of EPS.

Another promising application of the OPF has been the inclusion of switch and circuit breaker representation, hereinafter referred to as switchable branches, in power grid modeling. Through continuous functions, this representation has been used successfully in assessments of system reconfiguration [24] and in restoring the maximum loads in distribution systems [25]. The representation of these devices by integer variables is proposed in [26] to obtain the optimum dispatch and topology of an electric network to meet static loads. In some applications, the concept of switching transforms contingencies (criterion n − 1) into inequality constraints as a way of guaranteeing safety in the electrical system [27]. An OPF with a switching technique has also been used to alleviate transmission line congestion [28], in addition to improving the hosting of distributed generation systems [29] and minimizing losses [30] through network reconfiguration.

The conventional bus-branch representation of the power grid appears in the vast majority of OPF applications. However, developments related to the explicit representation of switches and circuit breakers in state estimation studies [31] have been incorporated into the power flow problem formulation [32], allowing for the processing of power grids at the substation level. In this type of modeling, the power flows through switchable branches are incorporated into the vectors of state variables together with the nodal complex voltages. This incorporation avoids possible numerical problems in the search for an optimal solution, motivated by the representation of impedances with very high or very low values. This is necessary to indicate the open and closed positions of these devices, respectively. As a consequence, this modeling allows for the direct determination of the power flow distribution through the components of substations, which is necessary for the analysis of electrical device loading, especially in circuit breakers. Preliminary studies demonstrated the feasibility of modeling at the substation level in the formulation of the OPF problem [33].

This work is inserted in this context and proposes an expeditious methodology capable of determining, for each hour of the day, the uncertainties in EPS electrical quantities owing to the uncertainties of wind and solar energy sources. The proposed approach consists of a Fuzzy Extended Deterministic OPF (FED-OPF) composed by two stages, which aims to explore the subjectivity of qualitative knowledge associated with these energy sources' behaviors. The first stage solves a deterministic OPF with extended formulation, referred to as ED-OPF, that is capable of processing modeled power grids at the substation level. In this first stage, the average values of the active power generated by wind and solar energy sources are used as input data together with the data referring to hourly load curves. In the second stage, the uncertainties of the EPS electric quantities are determined by the sum of the deterministic variables (determined in the first stage) with the fuzzy variables. The fuzzy variables are determined through qualitative sensitivity analysis applied to the postprocessing of the ED-OPF, considering maximum and minimum values for the active power of wind and solar energy sources. This procedure allows that the qualitative behavior of wind and solar energies sources be incorporated to the solution.

The proposed methodology allows for the assessment, in an expeditious way, of the influence of the uncertainties of the active power produced by wind and solar energy sources: (a) in the active power produced by conventional generators for energy reserve analysis, (b) in the nodal voltage magnitudes in order to comply with the limits imposed by regulatory authorities, and, mainly, (c) in the currents through circuit breakers. These results are essential to ascertain the correct sizing of the rated current of these devices.

#### **2. Fuzzy Extended Deterministic Optimal Power Flow**

Originating from an extension of the economic dispatch [34], the OPF allows for the determination of the optimum point of operation in relation to a predefined objective, while it complies with the operational limits imposed by power grids. Over the years, mainly owing to uncertainties in the generation and consumption of electric energy, the OPF was marked by a progression in the numerical techniques of solutions and in problem formulation. Among these techniques of solution, it is worth mentioning the Primal Dual Interior Points Method, initially proposed for the optimal dispatch of reactive power, which proved to be robust in solutions of large systems [35].

Based on these assumptions and considering the subjectivity of qualitative knowledge about the wind and solar energy sources behavior, this paper proposes a fuzzy extended deterministic optimal power flow (FED-OPF) capable of determining, for each hour of the day, the uncertainties of EPS electrical quantities owing to the uncertainties of these energy sources. These uncertainties are determined mainly for the currents through circuit breakers, to verify the correct sizing of these devices by the rated current, for the active power produced by conventional generators for energy reserve analysis and for the nodal voltage magnitudes in order to comply with the limits imposed by the regulatory authorities. The hourly analysis via the proposed FED-OPF is composed by two stages, as described in the following sections.

#### *2.1. Extended Deterministic Optimal Power Flow: First Stage*

In this work, the modeling of a power grid at the substation level is incorporated in the ED-OPF formulation. This modeling allows for the explicit representation of the switches and circuit breakers that make up the substation arrangements. This is different from bus-branch conventional modeling, which considers these arrangements as a single bus. To illustrate substation-level modeling, Figure 1 shows a simple five-bus system (a) where bus three started to rely on this modeling, resulting in an eleven-bus system (b).

**Figure 1.** (**a**) Bus-branch modeling and (**b**) substation level modeling.

With this new modeling, it is possible to obtain, in a direct way, information of electrical quantities through the substation components. This includes the power flows through circuit breakers, thus creating subsidies to confirm some technical specifications of these devices. As a consequence, this modeling excludes the need to use complementary analysis tools to obtain this information, as occurs with bus-branch conventional modeling.

It should be emphasized that the modeling at the substation level is applied only to previously selected substations, which are the substations through which it is desired to determine the power flow distribution and the influence of wind and solar sources. This means that only a small number, usually the substations connected to solar and wind sources, will be modeled at this level of detail.

In the following formulation, the details of the modeling of a power grid at the substation level and of the wind and solar power sources in the ED-OPF are presented:

$$\min \mathbf{C}(\mathbf{Pg}) = \mathbf{Pg}^{\dagger} \mathbf{Q} \mathbf{Pg} + \mathbf{b}^{\dagger} \mathbf{Pg} + \mathbf{co} \tag{1}$$

subject to:

$$\mathbf{P}\mathbf{g} + \overline{\mathbf{P}\mathbf{g}}\_{\text{win}} + \overline{\mathbf{P}\mathbf{g}}\_{\text{sol}} - \mathbf{P}\mathbf{d} = \text{real}\{\text{diag}(\mathbf{V}).(\mathbf{Y}\mathbf{V})^{\*}\} + \text{tfl} \tag{2}$$

$$\mathbf{Q}\mathbf{g} + \mathbf{B}\mathbf{s} - \mathbf{Q}\mathbf{d} = \text{imag}\{\text{diag}(\mathbf{V}).(\mathbf{Y}\mathbf{V})^{\*}\} + \text{uf}\mathbf{f} \tag{3}$$

$$\left(\boldsymbol{\Theta}\_{\bar{\mathbf{i}}\boldsymbol{\mathsf{j}}},\mathbf{t}\_{\bar{\mathbf{i}}}\right)^{\mathsf{P}\boldsymbol{\Theta}}=\mathbf{0}\tag{4}$$

$$\left(\mathbf{f}(\mathbf{V}\_{\vec{\text{ij}}}, \mathbf{u}\_{\vec{\text{ij}}})\right)^{\text{QV}} = 0 \tag{5}$$

$$\mathbf{P} \mathbf{g}\_{\min} \le \mathbf{P} \mathbf{g} \le \mathbf{P} \mathbf{g}\_{\max} \tag{6}$$

$$\mathbf{Q} \mathbf{g}\_{\min} \le \mathbf{Q} \mathbf{g} \le \mathbf{Q} \mathbf{g}\_{\max} \tag{7}$$

$$\mathbf{V}\_{\min} \le \mathbf{V} \le \mathbf{V}\_{\max} \tag{8}$$

$$\mathbf{a}\_{\min} \le \mathbf{a} \le \mathbf{a}\_{\max} \tag{9}$$

$$\mathbf{b}\_{\min} \le \mathbf{B} \mathbf{s} \le \mathbf{b}\_{\max} \tag{10}$$

$$
\mathfrak{sp}\_{\min} \stackrel{<}{\leq} \mathfrak{sp} \stackrel{<}{\leq} \mathfrak{sp}\_{\max} \tag{11}
$$

$$\mathbf{f}\mathbf{l}\_{\min} \le \mathbf{f}\mathbf{l} \le \mathbf{f}\mathbf{l}\_{\max} \tag{12}$$

Equation (1) represents the conventional objective function that minimizes the costs of active power generation **Pg** by conventional generators, where **Q** is the diagonal matrix with quadratic cost coefficients, **b** is the vector of linear cost coefficients and **co** is the constant cost vector. Equations (6) and (7) represent the minimum and maximum limit constraints of the active and reactive powers **Pg** and **Qg** produced by these generators, respectively. Equations (8) to (12) represent these limit constraints for the nodal voltage magnitudes **V**, taps of the transformers **a**, powers injected by static compensators **Bs**, angles of the shift transformers ϕ, and the active power flows on transmission lines **fl**, respectively.

The necessary adaptations in the ED-OPF formulation for the processing of power grids at the substation level are performed by extending the state vector so as to include the active (**t**) and reactive (**u**) power flow through each switchable branch along with the conventional nodal complex voltages. This procedure avoids the use of atypical values to represent the open and closed position of switchable branches in the problem formulation, eliminating the numeric problem such values would cause in the optimal solution search.

Equations (2) and (3) represents, respectively, the active and reactive bus power balance (or bus injected power), that is, the difference between the generated power and the demand at each bus. Please notice that, as a consequence of the current Kirchhoff law, the injected power at bus "k" is equal to the sum of the power flow through all adjacent branches of bus "k". When the substation level modeling is applied to the network, as proposed in this article, the sum of the power flow through adjacent branches must consider the set of adjacent conventional branches (transmission lines and transformers) as well as the set of adjacent switchable branches. Power flow through conventional branches are determined as usually, that is, as a function of the complex voltage. However, power flow through switchable branches are written directly as a function of the new state variables (**t** and **u**).

To incorporate the power flows through switchable branches in power balance equations (2) and (3) of the ED-OPF, the bus-switchable branches incidence matrix **A**bb is defined. With the number of rows equal to the number of buses of the system and the number of columns equal to the number of switchable branches, the element i-j of **A**bb is defined as


In (2), **t**fl represents the vector of the sum of the switchable branches active power flow in the buses. It can be written using the incidence matrix **A**bb as indicated in (13):

$$\text{tfl} = \mathbf{A}\_{\text{bb}} \mathbf{t} \tag{13}$$

Similarly, in (3), **u**fl represents the contribution, that is, the sum of the reactive power flows through switchable branches to the system buses, and can be written as:

$$\mathbf{u} \mathbf{f} = \mathbf{A}\_{\text{bb}} \mathbf{u} \tag{14}$$

where **t** and **u** are the new state variable vectors referring to the active and reactive power flows in switchable branches, respectively.

Equation (4) represents the active operational constraints that model the closed and open positions of the switches and circuit breakers. These constraints can be reformulated in function of specific incidence matrices for this purpose, as presented in (15):

$$\mathbf{f}(\boldsymbol{\Theta}, \mathbf{t})^{\mathrm{P}\boldsymbol{\Theta}} = \begin{bmatrix} \mathbf{A}\_{\mathrm{c}} & \mathbf{0} \\ \mathbf{0} & \mathbf{A}\_{\mathrm{o}} \end{bmatrix} \cdot \begin{bmatrix} \boldsymbol{\Theta} \\ \mathbf{t} \end{bmatrix} = \begin{bmatrix} \mathbf{A}\_{\mathrm{c}} \boldsymbol{\Theta} \\ \mathbf{A}\_{\mathrm{o}} \mathbf{t} \end{bmatrix} \tag{15}$$

The incidence matrix **A**<sup>c</sup> has a number of lines equal to the number of closed circuit breakers and a number of columns equal to the number of total buses of the system. This matrix represents the operational constraints of closed circuit breakers of the type θ<sup>i</sup> − θ<sup>j</sup> = 0, where θ is the angle of the complex voltages. Thus, the values of this matrix are equal to 1 in column "i" and equal to -1 in column "j" of the line corresponding to the closed circuit breaker i-j.

The incidence matrix **A**o has a number of lines equal to the number of open circuit breakers and a number of columns equal to the total number of circuit breakers. This matrix represents the operational constraints of open circuit breakers of type **t**ij = 0. The values of this matrix will only be different from zero and are necessarily equal to 1 in the line corresponding to open circuit breaker i-j and the column corresponding to the state variable (active power flow) associated with this circuit breaker.

A similar reformulation can be performed on the reactive operating restrictions of (5), as shown in (16):

$$\mathbf{f}(\mathbf{V}, \mathbf{u})^{\mathrm{QV}} = \begin{bmatrix} \mathbf{A}\_{\mathrm{c}} & \mathbf{0} \\ \mathbf{0} & \mathbf{A}\_{\mathrm{0}} \end{bmatrix} \cdot \begin{bmatrix} \mathbf{V} \\ \mathbf{u} \end{bmatrix} = \begin{bmatrix} \mathbf{A}\_{\mathrm{c}} \mathbf{V} \\ \mathbf{A}\_{\mathrm{0}} \mathbf{u} \end{bmatrix} \tag{16}$$

From (16), it is possible to obtain the operational restrictions of the closed circuit breakers of type Vi − Vj = 0, in addition to the operational restrictions of the open circuit breakers of type **u**ij = 0.

Regarding the representation of alternative sources of energy, (2) considers in the balance of the active power the average values of the active power produced by wind and solar energy sources **Pg**win and **Pg**sol, respectively. The average values can be determined from real historic data, as discussed in the Section 3.

The solution of the proposed extended deterministic optimal power flow is achieved by the Primal Dual Interior Points Method [35] as associated with the Lagrange technique, the conditions of Karush Kuhn Tucker, and the Newton Method.

This solution defines the deterministic variables of EPS electrical quantities as function of average values of the active power produced by wind and solar energy sources and demands the traditional expensive computational time of solutions involving OPFs.

#### *2.2. Electrical Quantity Uncertainties of EPS: Second Stage*

The information regarding the uncertainties of generation and consumption of electric energy in the OPF is usually represented by probabilistic approaches through probability density functions [4,36,37]. However, these approaches based on quantitative information are not capable of representing the qualitative information that resides in human knowledge. With the fuzzy technique, the incomplete, vague, and qualitative knowledge present in human beings began to be represented by mathematical models in several areas of knowledge, including in OPFs [15].

In addition to not considering qualitative information, usually the processing of probabilistic data in the OPF demands a high computational effort. However, in several practical situations, expeditious and even approximate responses of the electric quantities of the power systems are enough for the analysts of these systems. An example of such situations is the investigation of the necessity of replacing circuit breakers with equipment that has a greater rated current capacity, usually common in expansions and changes in EPS.

Based on the above discussion, this work proposes an hourly evaluation of the electrical quantities of the electric system in relation to the intermittent behavior of wind and solar energy sources, based on the association illustrated in the Figure 2. In this illustration the area delimited by the mean and +/−3 standard deviations of the Probability Density Function (PDF) normal distribution of the active power produced by wind and solar energy sources is associated with the area delimited by the membership function of a fuzzy triangular model. This approach allows composing, through fuzzy variables and a real brief historical data about these energy sources, an hourly qualitative knowledge about the average, minimum, and maximum values of powers produced by these sources. Section 3 exemplifies the composition of this qualitative knowledge.

**Figure 2.** Probability Density Function (PDF) and fuzzy triangular model for wind and solar energy sources.

The Figure 2 shows the Probability Density Function (PDF) (in black), where: *f*(**Pg**win): probability density function of the active power produced by wind energy sources; **Pg**win: hourly average values of the active power produced by wind energy sources; *σ*: standard deviation of the active power produced by wind energy sources.

For the fuzzy triangular model (in blue), the Figure 2 shows: μ(**Pg**win): membership function of the active power produced by wind energy sources; Δ**Pg**win: hourly variations in relation to the average values of the active power produced by wind energy sources, here called fuzzy variables. For the solar energy sources, we can use the same definition, changing the subscript "win" to "sol".

As the qualitative average values of the active power produced by these sources **Pg**win and **Pg**sol were considered in the optimal solution of the ED-OPF of the first stage (Equation (2)), it is possible to define the fuzzy variables of the other system electrical quantities as a function of the fuzzy variables Δ**Pg**win. and Δ**Pg**sol. For this purpose, let us consider the Lagrangian function of the problem formulation (1)–(12). The first derivative (or gradient) of this function in the optimal solution must be null, in order to comply with the optimality conditions of Karush Kuhn Tucker. This gradient is represented by (17):

$$
\rho \left( \mathbf{P} \mathbf{g}\_{\prime} \, \overline{\mathbf{P}} \overline{\mathbf{g}}\_{\text{win}} , \, \overline{\mathbf{P}} \overline{\mathbf{g}}\_{\text{sol}} , \mathbf{V} , \mathbf{t} , \mathbf{y} \right) = 0 \tag{17}
$$

where, **Pg** is the vector of the active power produced by conventional generators; **Pg**win and **Pg**sol are the vectors of the qualitative average values of the active power produced by wind and solar energy sources, respectively; **V** is the vector of the nodal voltage magnitudes; **t** is the vector of the active power flows through circuit breakers; and **y** is the vector of the other electrical quantities of the system, including the dual variables.

Based on this, the expansion of (17) in Taylor series, in the directions of fuzzy variables Δ**Pg**win and Δ**Pg**sol, still in the optimal solution, can be formulated by:

$$\begin{split} \rho \left( \mathbf{P} \mathbf{g} + \boldsymbol{\Lambda} \mathbf{P} \mathbf{g}, \, \overline{\mathbf{P} \mathbf{g}}\_{\text{virin}} + \overline{\boldsymbol{\Lambda} \mathbf{P} \mathbf{g}}\_{\text{virin}}, \, \overline{\mathbf{P} \mathbf{g}}\_{\text{sol}} + \overline{\boldsymbol{\Lambda} \mathbf{P} \mathbf{g}}\_{\text{sol}}, \, \mathbf{V} + \boldsymbol{\Lambda} \mathbf{V}, \mathbf{t} + \boldsymbol{\Lambda} \mathbf{t}, \, \mathbf{y} + \boldsymbol{\Lambda} \mathbf{y} \right) \\ = \rho \left( \mathbf{P} \mathbf{g}, \, \overline{\mathbf{P} \mathbf{g}}\_{\text{vair}}, \overline{\mathbf{P} \mathbf{g}}\_{\text{sol}}, \mathbf{V}, \mathbf{t}, \mathbf{y} \right) + \frac{\partial \mathbf{p}}{\partial \mathbf{P} \mathbf{g}} \boldsymbol{\Lambda} \mathbf{P} \mathbf{g} + \frac{\partial \mathbf{p}}{\partial \mathbf{P} \mathbf{g}}\_{\text{sol}} \boldsymbol{\Lambda} \overline{\mathbf{P} \mathbf{g}}\_{\text{vair}} + \frac{\partial \mathbf{p}}{\partial \mathbf{P} \mathbf{g}}\_{\text{sol}} \boldsymbol{\Lambda} \overline{\mathbf{P} \mathbf{g}}\_{\text{sol}} + \frac{\partial \mathbf{p}}{\partial \mathbf{V}} \boldsymbol{\Lambda} \mathbf{V} + \frac{\partial \mathbf{p}}{\partial \mathbf{I}} \boldsymbol{\Lambda} \mathbf{t} + \frac{\partial \mathbf{p}}{\partial \mathbf{V}} \boldsymbol{\Lambda} \mathbf{y} \end{split} \tag{18}$$

This relationship implies that:

$$
\begin{bmatrix}
\frac{\partial \mathbf{p}}{\partial \mathbf{P} \mathbf{g}} \ \frac{\partial \mathbf{p}}{\partial \mathbf{V}} \frac{\partial \mathbf{p}}{\partial \mathbf{t}} \ \frac{\partial \mathbf{p}}{\partial \mathbf{y}}
\end{bmatrix}
\begin{bmatrix}
\Delta \mathbf{P} \mathbf{g} \\
\Delta \mathbf{V} \\
\Delta \mathbf{t} \\
\Delta \mathbf{y}
\end{bmatrix} = -\frac{\partial \mathbf{p}}{\partial \overline{\mathbf{P} \mathbf{g}}\_{\text{win}}} \Delta \overline{\mathbf{P} \mathbf{g}}\_{\text{win}} - \frac{\partial \mathbf{p}}{\partial \overline{\mathbf{P} \mathbf{g}}\_{\text{sol}}} \Delta \overline{\mathbf{P} \mathbf{g}}\_{\text{sol}}\tag{19}
$$

Or that:

$$
\Delta \mathbf{W}' \begin{bmatrix} \Delta \mathbf{P} \mathbf{g} \\ \Delta \mathbf{V} \\ \Delta \mathbf{t} \\ \Delta \mathbf{y} \end{bmatrix} = -\frac{\partial \mathbf{p}}{\partial \overline{\mathbf{P}} \mathbf{g}\_{\text{win}}} \Delta \overline{\mathbf{P}} \overline{\mathbf{g}}\_{\text{win}} - \frac{\partial \mathbf{p}}{\partial \overline{\mathbf{P}} \mathbf{g}\_{\text{sol}}} \Delta \overline{\mathbf{P}} \overline{\mathbf{g}}\_{\text{sol}} \tag{20}
$$

Defining **W** as a Hessian matrix, the equation can be rewritten as:

$$
\begin{bmatrix}
\Delta \mathbf{P} \mathbf{g} \\
\Delta \mathbf{V} \\
\Delta \mathbf{t} \\
\Delta \mathbf{y}
\end{bmatrix} = \mathcal{W}^{-1} \left( -\frac{\partial \mathbf{p}}{\partial \overline{\mathbf{P}} \mathbf{g}\_{\text{win}}} \Delta \overline{\mathbf{P}} \overline{\mathbf{g}}\_{\text{win}} - \frac{\partial \mathbf{p}}{\partial \overline{\mathbf{P}} \mathbf{g}\_{\text{col}}} \Delta \overline{\mathbf{P}} \overline{\mathbf{g}}\_{\text{col}} \right) \tag{21}
$$

Equation (21) defines the fuzzy variables of the active power produced by conventional generators Δ**Pg**, of the nodal voltage magnitudes Δ**V**, of the active power flows through switches and circuit breakers Δ**t**, and of the other system quantities Δ**y**. With these fuzzy variables, it is possible to determine the uncertainties of the electrical quantities of the EPS. These uncertainties result from the sum of the deterministic variables calculated in the optimal ED-OPF solution as determined in the first stage, with the fuzzy variables determined in this stage, as shown below:

$$\mathbf{P} \mathbf{g}\_{\text{unccertain}} = \mathbf{P} \mathbf{g} + \ \Lambda \mathbf{P} \mathbf{g} \tag{22}$$

$$\mathbf{V}\_{\text{unccertain}} = \mathbf{V} + \Delta \mathbf{V} \tag{23}$$

**t**uncertain = **t** + Δ**t** (24)

$$\mathbf{y}\_{\text{uncertain}} = \mathbf{y} + \Delta \mathbf{y} \tag{25}$$

where **Pg**uncertain are the uncertainties regarding the active power produced by conventional generators, necessary for energy reserve analysis; **V**uncertain are the uncertainties regarding the nodal voltage magnitudes, necessary to analyze the accomplishment of the limits imposed by the regulatory authorities; **t**uncertain are the uncertainties regarding the active power flows through circuit breakers, necessary to verify the correct sizing of these devices by the rated current; and **y**uncertain are the uncertainties regarding the other electrical quantities that can be used according to necessity.

In summary, the uncertainties of EPS electrical quantities are determined, in the second stage, by the sum of deterministic variables (determined in the first stage) and fuzzy variables. The fuzzy variables are determined through qualitative sensitivity analysis (Equations (17)–(21)) applied to the postprocessing of the ED-OPF optimal solution, as function of variation (maximum and minimum values) of the active power of wind and solar energy sources, always considering an hourly based analysis.

The fuzzy variables participate only in the postprocessing stage of the OPF and avoid the incorporation of these uncertainties into the iterative process of the OPF solution. The result is a significant reduction in the computational time when compared with the first stage. The second stage provides rapid and authentic responses (as discussed in next section) for a general analysis of the electrical quantities of the system as a function of the variation of wind and solar energy sources.

#### **3. Test System Data**

To implement the proposed methodology, a real power network composed of 139 buses that comprises the northeast Brazilian coastal region was used. The system contains buses at different voltage levels ranging from 13.8 kV to 500 kV. In addition, the available data includes real hourly load curves, static compensators, capacitor banks, bus and line reactors, previously selected substations modeled at the substation level, and large hydroelectric power plants: Sobradinho (1.00 GW), Xingó (3.6 GW), Luiz Gonzaga (1.5 GW), and Paulo Afonso Complex (2.4 GW). The power system has approximately 400 MW of thermoelectric power plants.

Additionally, the power system has 19 wind farms that represent approximately 2 GW of installed power. Figure 3 illustrates the daily variation of active power produced by these wind farms in the period between June 2016 and July 2017. All of this information corresponds to daily real data properly processed and grouped to subsidize the analysis of the proposed tool in the most realistic way possible.

**Figure 3.** Maximum, average, and minimum daily values of active power (MW) wind farms.

To illustrate the proposed qualitative definition of maximum, average, and minimum of active power of the renewable energy sources, the Brisa Potiguar wind farm (BP) were chosen. Figure 4 presents the active power average values **Pg**win of BP for each hour of the day, which was determined from the available real historic data.

**Figure 4.** Brisa Potiguar wind farm daily behavior.

These values were used in the first stage (Section 2.1—Equation (2)) of the proposal methodology, in order to define the deterministic variables of EPS electrical quantities.

According to real historical data, the maximum and minimum values of active power of BP are defined for each hour of the day. These values correspond to the variation of +/−30% of the average values of these powers, which are represented the fuzzy variables Δ**Pg**win(see Figure 2). These values are used in the second stage (Section 2.2—Equations (17)–(21)) of the proposal methodology, in order to determine the fuzzy variables of EPS electrical quantities. The uncertainties of EPS electrical quantities are then computed by the sum of deterministic variables, determined in the first stage, and fuzzy variables, determined in the second stage, as discussed in Section 2.2—Equations (22–25).

In summary, the average values of active power of wind and solar energy sources are determined through the brief real historic data. The maximum and minimum values of these active powers are defined by the hourly variation of +/−30% of the average value. It should be emphasized that the whole proposed strategy assumes an hourly based analysis, such that fuzzy triangular model (based in the normal PDF) is used to model the hourly variation of the wind and solar generators in the second stage of the proposed methodology.

The same methodology is applied to active power produced by solar energy sources, which are randomly distributed in six of the 69 kV buses. The solar energy participates in the formulation of the problem only in the period of solar incidence between 6–18 h, and the average value of these powers **Pg**sol is dimensioned in 10% of the active power demanded in these buses. The fuzzy variables Δ**Pg**sol, related to the active power produced by solar energy sources are also defined by values equivalent to +/−30% of the average values.

#### **4. Simulations and Results**

The developments presented in this work were carried out using the MATLAB computational tool. The first and most important result consists of authenticating the proposed methodology. In order to validate the proposed approach, exhaustive solutions of conventional OPF considering each increase/decrease variation for wind and solar sources, were used.

Figure 5 illustrates this validation for the active power of the hydroelectric plant of Xingó.

**Figure 5.** Authentication—active power of Xingó hydroelectric power plant.

The "Original" trace of Figure 5 represents the deterministic values of the active power of the hydroelectric plant of Xingó, obtained with the ED-OPF solution, as described in Section 2.1. In this case, the qualitative average values of the active power of the wind and solar energy sources were used in the power balance equations. The "Wind and Solar Increment" trace represents the uncertainties of the active power of the hydroelectric plant of Xingó. These uncertainties are determined by the sum of determinist values ("Original" trace) with fuzzy values. In this case, according to the technique described in Section 2.2, the fuzzy values were determined through qualitative sensitivity analysis applied to the postprocessing of the ED-OPF, as function of 30% (maximum values) increase in the active power of wind and solar energy sources. The "Wind and Solar Decrement" trace also represents the uncertainties of the active power of the hydroelectric plant of Xingó, but now considering a decrease of 30% (minimum values) in the active power of these energy sources.

Simulations were performed in a conventional OPF to authenticate the proposed methodology. The "Original Validation" trace represents the values obtained in the solution of a conventional OPF that considers in the power balance equations the qualitative average values of the active power of the wind and solar energy sources. The "Wind and Solar Increment Validation" trace represents the values obtained in the conventional OPF that considers in the active power balance equations the qualitative average values of the active power of these energy sources increased by 30%, while the "Wind and Solar Decrement Validation" trace represents the values obtained in the conventional OPF when in these equations the qualitative average values of the active power of these energy sources is reduced in 30%. The adherence of the results indicates the effectiveness of the proposed methodology in obtaining the desired results in an expeditious way, as proposed.

This authentication is also successfully achieved for the other electrical quantities of the system. The following sections present the definition of the uncertainties of the power flows through circuit breakers, of the active power produced by conventional generators, and of the nodal voltage magnitudes as a function of the uncertainties of the active power produced by wind and solar energy sources. These analyses are academic and have no intention of interfering with the operation or expansion of the actual system.

Similar proposals, that reach solution with reduced computational time as proposed in our paper, and that simultaneously contemplates in the OPF, renewable energies, switching (substation model level), and qualitative information, were not found in the literature. These characteristics make us confident that the proposed approach presents innovative ideas and presents itself as a promising tool for analyzing the impact of intermittent energy sources in the electrical quantities of the power systems.

#### *4.1. Power Flows through Circuit Breakers*

The verification of the correct sizing of the rated current of circuit breakers, often motivated by expansions or alterations in power systems, has attracted the attention of analysts of these systems

in different parts of the world. In these analyses, the representation of the substation of interest at the substation level, as proposed in this paper, allows for the assessment in a direct way of the power flows through these devices. As most analytical tools consider bus-branch modeling, this analysis requires the adoption of additional procedures that overload the area professionals and hamper this task. The incorporation of network modeling at the substation level in the OPF eliminates this step, thus facilitating and more effectively subsidizing the analysis work by these professionals.

To illustrate the effectiveness of this attribute of the proposed methodology, the Sobradinho 230 kV substation was chosen for modeling at the substation level owing to its connections with wind and solar power sources, as well as sheltering an important hydroelectric plant of the same name with 1.05 GW of installed power. Figure 6 illustrates the configuration of this substation, including the most common directions of power flows through the circuit breakers of this installation.

**Figure 6.** Power flow distribution at Sobradinho 230 kV substation.

Table 1 presents the most relevant results related to this analysis involving the Sobradinho 230 kV substation. Column 1 indicates the hour at which the variation was assessed (corresponding to the solar incidence period), while the second to fifth columns show the variations in the active power flows through the indicated circuit breaker when the active power produced by wind and solar sources varied from +/−30% (always in relation to the average values according to Section 2.1).


**Table 1.** Active Power Flow Variations (MW) Through Circuit Breakers of Sobradinho 230 kV Substation.

The results indicate that the active power flows through the circuit breakers of this substation were influenced by the simultaneous generation of wind and solar energy sources. The exception was with the active power flows through circuit breaker 159, which were influenced only by the wind energy source since this device was connected to the Pedra Branca 230 kV wind farm. Similar analyses can be performed for the active power flows through circuit breakers 160 and 163, which connect the solar source, and therefore were influenced only by this source of energy.

Once the power flows through the circuit breakers of the Sobradinho 230 kV substation were determined, the maximum currents through these devices were calculated as a function of the apparent powers. From this calculation, it was observed that the maximum ratio between these currents and the rated current was 22.2% for circuit breaker 159, ensuring the correct dimensioning of these devices by the rated current.

The João Câmara II 138 kV bus was also modeled at the substation level owing to the connection of approximately 700 MW of installed power by wind farms. Figure 7 shows the bus configuration and indicates the most common orientation of power flows through the circuit breakers of this substation.

**Figure 7.** Power flow distribution at João Câmara II 138 kV substation.

Table 2 lists the maximum current values through the main circuit breakers of this substation, as well as the proportion of these values in relation to their rated values. These analyses are performed for the original condition, that is, without the variation of the average values of the active power of the wind sources, and with a variation of +/−30% in relation to these values.


**Table 2.** Maximum Currents (A) Through Circuit Breakers of João Câmara II 138 kV Substation.

The results indicate that the effect of the wind power sources was very significant, taking the maximum current values through these devices near the rated values. For example, the maximum current through circuit breaker 194 reached 71.6% of its rated current value when exposed to the uncertainties of the active power produced by the wind farms connected to the bus under study. The effect of the solar energy sources in this substation was negligible.

It should be noted that this case considers a normal operating condition. Logically, under conditions of contingency (criterion n − 1) or expansions in the system considered, the currents passing through these circuit breakers could be closer to the nominal values and may even exceed them. This type of analysis can be carried out by the proposed methodology and allows for the assessment of the necessity of replacing these devices with equipment that has a greater rated current capacity.

#### *4.2. Active Power: Conventional Generators*

Another relevant issue associated with the uncertainties of the active power generated by wind and solar energy sources is related to the analysis of their impact on the established powers for conventional generators. Although it is not advisable to ignore the importance of intermittent sources in reducing the rotating reserve, a careful analysis is also needed to avoid compromising the energy supply to consumers in conditions of unexpected reductions of the share of these energy sources in the energy matrix.

Table 3 lists, for some hours of the day, the variation in megawatts from the average values of the active power produced by the Paulo Afonso IV hydroelectric power plant. These values were obtained as a function of the variation of +/−30% from the average values of the active power generated by wind and solar energy sources.


**Table 3.** Variations (MW) at the Paulo Afonso IV Hydroelectric Power Plant.

These and other results obtained (omitted because of lack of space) show that all conventional generators of the system were sensitized by both the solar and wind energy uncertainties except for the generators of the thermoelectric plants that were not dispatched owing to high cost.

In this context, Table 4 shows the influence of the uncertainties of wind and solar energy on the energy reserve of these conventional generators.


**Table 4.** Available Powers: Conventional Generators.

The results presented show a slight reduction in the available powers of generators, which was necessary to increase energy production to supply the reduction in the share of wind and solar sources. The Luiz Gonzaga hydroelectric power plant reached its active power limit with a reduction of 30% in the share of these energy sources. The system under analysis contemplates an installed capacity of approximately 20 GW distributed between the hydroelectric and thermoelectric plants and equivalent

generators, in addition to the hourly demand of around 7 GW. Considering energy generation planning based on an operational reserve equivalent to 5% of the demand, the variations of +/−30% of the active power of the wind and solar energy sources guarantee an operational reserve of 350 MW.

#### *4.3. Nodal Voltage Magnitudes*

Voltage control is one of the main concerns when there are expansions and changes in energy systems, as it must comply with the limits imposed by regulatory authorities. Thus, for the test system of this work, the uncertainties of the nodal voltage magnitudes were determined as a function of the uncertainty of +/−30% in relation to the average value of the active power of the wind and solar energy sources. Figure 8 illustrates the behavior throughout the day of the nodal voltage magnitudes in the 500 kV Pau de Ferro substation, an important transmission trunk of the system, owing to the combination of the uncertainties of wind and solar energy sources.

**Figure 8.** Uncertainties in nodal voltage magnitudes of 500 kV bus of Pau de Ferro substation.

The values obtained for the nodal voltage magnitudes of all buses of the test system are within the limits established by the Brazilian regulatory authority, according to Table 5 [38]. These results show once again the relevance of the proposed methodology that allows for the inference of questions of this order without the need to carry out studies and exhaustive simulations involving the intermittent sources of energy.


**Table 5.** Voltage magnitude admissible limits in the Brazilian electrical system.

#### *4.4. Computational Time*

Table 6 presents the computational times demanded in the execution of first and second stages of the proposed methodology.


**Table 6.** Computational time (seconds).

According to the results presented in Table 6 it is possible to verify the significant computational time reduction in the second stage.

Based on the subjectivity of the qualitative knowledge about the wind and solar energy sources behavior, the second stage of the proposal methodology provides for the analysts fast answers and with good authenticity about electrical quantities of these systems, as indicated in Figure 5.

#### **5. Conclusions**

The combination of fuzzy logic with a Taylor series in the postprocessing of the Extended Deterministic Optimum Power Flow allows for the assessment of the impact of the active power uncertainties produced by the wind and solar energy sources in the electrical quantities of energy systems, for each hour of the day. In an expeditious way, based in a subjective investigation that focuses on the qualitative character of these energy sources' behaviors, the proposed methodology allows for the assessment of the impact of these uncertainties on the active power produced by conventional generators. This makes it possible to conduct energy reserve analyses. These impacts are also assessed on the nodal voltage magnitudes. Thus, it is possible to verify if the limits reached comply with the limits of the voltages imposed by the regulatory authorities.

In addition, a representation of the power grid at the substation level made it possible to identify the impact of the variations of the wind and solar energy sources on the distribution of power flows and consequently of the currents through the circuit breakers of substations. This model allows for the evaluation of the technical specifications of such devices, such as rated currents. In face of this new reality of intermittent energy generation systems, this topic has attracted the attention of systems analysts from various parts of the world owing to the need for expansions and changes in EPS.

The results presented highlight the importance of energy reserve analyses and of correct technical specifications of circuit breakers in addition of voltages control, as a function of the forecasted growth of wind and solar energy sources in energy matrices. These results reinforce the relevance of analytical tools that provide professionals the ability to perform these tasks, in an expeditious way, as in the case of the methodology proposed in this work. Future studies are under investigation to assess the impact of such approach in contingencies analysis, including wind farm groups contingencies.

#### **6. Article Contribution**

The main contribution of this article is to offer power systems analysts a tool that allows the rapid assessment of the hourly behavior of the electrical quantities of these systems as a function of the variation of wind and solar energy sources. It is not a tool that involves data accuracy, but rather a tool that provides enough and rapid responses to decision making by those analysts. The proposed tool is original and suitable mainly for analysis of currents through circuit breakers, a theme that has attracted the attention of analysts from various parts of the world, due to the need for expansions and topology changes in power systems. The variations in the power injected by these energy sources can increase the value of the currents through the circuit breakers, resulting in an alert to the analysts about the need for replacement with equipment with higher nominal current value.

The assessment of the impact of the variation of the wind and solar energy sources not only on the currents through the circuit breakers, but also on the reserve energy and the nodal voltage magnitudes is reached through the solution of fuzzy extended deterministic optimal power flow, composed of two stages.

As discussed before, in the first stage, the extended deterministic OPF (ED-OPF) solution provides the deterministic variables of the power system electrical quantities. For this stage, we would like to emphasize that the formulation of ED-OPF contemplates an original modeling at the substation level, which enables the direct assessment of the currents through the circuit breakers. In addition, this formulation contemplates the hourly average values of the active powers of wind and solar energy sources in the active power balance equations. The first stage solution demands the traditional expensive computational time of solutions involving OPFs.

In the second stage, the uncertainties of the power system electrical quantities are obtained considering an hourly based analysis. These uncertainties are determined by the sum of determinist variables (first stage) with fuzzy variables. The fuzzy variables are obtained when the maximum and minimum values of these energy sources (corresponding to +/−30 % of average values of these energy sources) are applied to the ED-OPF solution, through a qualitative sensitivity analysis. The second stage presents lower computational time when compared to the first stage.

With the proposed tool, power system analysts can run ED-OPF once, which demands longer computational time and work more frequently with the second stage, which demands lower computational time. The second stage provides rapid and authentic responses for a general analysis of the electrical quantities of the system as a function of the variation of wind and solar energy sources.

**Author Contributions:** All the authors contributed to the realization, analysis of the data, as well as the writing of the paper.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are grateful for the support of the Brazilian electric utilities Companhia Paranaense de Energia (COPEL) (Paraná, Brazil) and Companhia Hidroelétrica do São Francisco (CHESF) (Pernambuco, Brazil).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Optimal Placement of UHF Sensors for Accurate Localization of Partial Discharge Source in GIS**

**Rui Liang 1,2,\*, Shenglei Wu 3, Peng Chi 1, Nan Peng <sup>1</sup> and Yi Li <sup>1</sup>**


Received: 12 February 2019; Accepted: 21 March 2019; Published: 26 March 2019

**Abstract:** This paper proposes an optimal placement model of ultra-high frequency (UHF) sensors for accurate location of partial discharge (PD) in gas-insulated switchgear (GIS). The model is based on 0-1 program in consideration of the attenuation influence on the propagation of electromagnetic (EM) waves generated by PD in GIS. the optimal placement plan improves the economy, observability, and accuracy of PD locating. After synchronously acquiring the time of the initial EM waves reaching each UHF sensor, PD occurring time can be obtained. Then, initial locating results can be acquired by using the Euclidean distance measuring method and the extended time difference of arriving (TDOA) location method. With the information of all UHF sensors and the inherent topological structure of GIS, the locating accuracy can be further improved. The method is verified by experiment, showing that the method can avoid the influence of false information and obtain higher locating accuracy by revising initial locating results.

**Keywords:** gas insulated switchgear; partial discharges; electromagnetic wave propagation; optimization; UHF measurement

#### **1. Introduction**

With the wide use of gas-insulated switchgear (GIS), various failures occur with the growth of service time, and the insulation defects in GIS may lead to partial discharge (PD). PD location can be realized via ultrasonic waves and ultra-high frequency (UHF) electromagnetic (EM) waves [1–3]. The UHF method is an efficient PD detection technology that detects the EM wave generated by PD based on UHF sensors [4–8]. The UHF method has the advantages of high sensitivity and high signal-to-noise ratio. The position of PD can be located by analyzing and processing the EM signal received by the sensors [9–11]. UHF location methods are mainly based on time difference of arriving, time of arriving, angle of arriving, and received signal strength [12,13].

The spatial placement of UHF sensors is the key to PD locating. Studies have been undertaken to improve sensor economy and observability. Though some exploratory research works have been carried out in these fields, most of them focused on transformer or conventional switchgear [14–18]. For instance, a spatial array of UHF sensors is used for transformer PD location. However, a transformer is a multisurface system of which the internal structure is totally different from GIS. The spatial array aims to overcome the reflection and diffraction of electromagnetic waves. Therefore, the methods for transformers and switchgear are not suitable for GIS. The topological structures of GIS and the attenuation of the EM wave especially propagating through the insulators and L-type and T-type structures should be considered more. A few studies have investigated the topological structures of GIS but these do not consider the optimal placement of sensors [19–23]. The arriving time difference

of EM waves was used to locate the PD in GIS based on UHF sensors, but its accuracy should be improved [24–26]. This paper is devoted to solving these problems in existing research.

In this paper, the GIS space structure is abstracted into an undirected graph and a 0-1 programming model is established based on the optimization of the optimal economy and the maximum observability. Then, an arrangement plan of UHF sensors can be obtained. With the optimal placement plan, the accurate location method of PD based on extended time difference of arriving (TDOA) is studied. The experimental result indicates that the proposed method can improve the location accuracy and possesses fine fault tolerance.

#### **2. Principle of Method**

#### *2.1. Partial Discharge in GIS*

GIS, which includes a circuit breaker, disconnector, grounding switch, transformer, and bus, has been widely used due to its low maintenance, high reliability, and compact size. PD refers to the discharge phenomenon occurring in the partial area of the insulator which does not penetrate the conductor under voltage. Under a strong electric field, the insulation defects of GIS will lead to PD, accompanied by the generation of electricity, light, heat, sound, ozone, and nitric oxide, which will corrode the insulating materials. In addition, the charged particles produced by PD will impact the insulating materials. PD is a sign of insulation degradation in GIS, which endangers the safety of equipment and even the power system.

The breakdown field strength of insulators in power equipment is very high. When PD occurs in a small range, the breakdown process is very fast, which produces steep pulse current. Its rising time is less than 1 ns, and ultra-high frequency (UHF) electromagnetic waves are excited. The principle of the UHF method for PD detection is to detect the UHF electromagnetic waves (300 MHz-3 GHz) generated by PD in power equipment by UHF sensor. Then, the information of PD can be obtained. Built-in UHF sensors and external UHF sensors are usually used for UHF methods.

Because corona interference is mainly concentrated below the 300 MHz frequency band, the UHF method can effectively avoid corona interference. It has high sensitivity and anti-interference ability. Naturally, the optimal placement of UHF sensors and location method become the key to partial discharge locating.

#### *2.2. Optimal Placement of UHF Sensors*

The optimal placement of UHF sensors should satisfy the economy of scale and provide maximum observability. A necessary requirement for successful locating is that at least two sensors detect an effective discharge EM wave, and the EM wave must pass through both ends of the discharge section. If the EM wave only passes through the head or the end of the detecting section, the locating will fail. 1 represents that the sensor needs to be installed and 0 represents the sensor doesn't need to be installed. Therefore, the placement of UHF sensors can be abstracted as a 0-1 programming model.

The EM wave will be attenuated when passing through the insulators and L-type and T-type structures in GIS. The vertical branches of L-type structures and T-type structures significantly attenuate the EM wave, and the attenuation caused by the horizontal branches of T-type structures and insulators comes after. In order to improve the detection sensitivity and locating reliability, the L-type structures and the T-type structures in GIS should be taken into consideration.

In topology, the number of branches connected to the node is called the degree (*d*) of the node. Obviously, the degree of an L-type node (node 4 in Figure 1) is 2 and the degree of a T-type node (node 4 in Figure 2) is 3. For general GIS, the node with the largest degree is the cross node, of which the degree is 4, and the node with the smallest degree is the terminal node, of which the degree is 1. The cross node can be considered as the combination of two L-type nodes. To effectively detect the EM wave of PD, the EM wave should not pass through more than one L-type node, one cross node, or one

T-type node. These three nodes are collectively called the nonterminal node, of which the degree is greater than or equal to 2.

**Figure 1.** L-type structure of gas-insulated switchgear (GIS).

**Figure 2.** T-type structure of GIS.

The number of branches included in the shortest path of node *i* to node *j* is defined as the shortest path length *Pij*. Under the condition of installing sensors at the terminal node (subject to *xt* = 1), a suitable nonterminal node should be selected to install a sensor so that the EM wave passes through, at most, one nonterminal node before reaching the sensors. Therefore, the maximum number of nonconfigurable nodes between any adjacent configurable nodes is 1, and *Pij* should be no more than 2 (subject to *Pij* = 2). An optimization model can be obtained:

$$\begin{array}{l}\min \sum\_{k=1}^{N} x\_k\\ \text{s.t. } P\_{ij} \le 2\\ x\_t = 1 \end{array} \tag{1}$$

where *t* is terminal node, *i* and *j* are adjacent nodes, and *xk* is the node in which the sensor should be installed. Because the configurable node is unknown, the inequality constraint of the shortest path which contains variable *xk* cannot be obtained directly.

The condition above can be described as follows: at least one sensor should be installed on those nodes of which the shortest path length to node *k* is 1 or 2.

For a node *<sup>k</sup>*, an adjacent node matrix *<sup>V</sup><sup>k</sup>* is proposed. *<sup>V</sup><sup>k</sup>* is *<sup>M</sup>* × *<sup>N</sup>* matrix, where *<sup>N</sup>* is the number of nodes and *M* is the number of the nodes, of which the shortest path length to the node *k* is 2. The nodes of which the shortest path length is 1 don't need to be considered because the shortest path from these *<sup>M</sup>* nodes to the node *<sup>k</sup>* includes them. For a row vector *<sup>v</sup>* of *<sup>V</sup><sup>k</sup>* and *<sup>j</sup>* (1 ≤ *<sup>j</sup>* ≤ *<sup>N</sup>*), when *Pkj* = 1 or *Pkj* = 2, *vj* = 1 and other elements are 0. Therefore, it forms the constraint that at least one sensor should be installed on the path from the node *k* to the node of which the shortest path length is 2. *<sup>V</sup>*1~*<sup>N</sup>* is arranged in columns to form a new matrix *<sup>V</sup>M*×*N*, and it can be found as follows:

$$\mathbf{V}\_{\mathbf{M}\times\mathbf{N}}\mathbf{X}\_{\mathbf{N}\times\mathbf{1}}\geq\mathbf{I}\_{\mathbf{M}\times\mathbf{1}}\tag{2}$$

where *<sup>X</sup>*<sup>T</sup> = [*x*1, *<sup>x</sup>*2,... , *xN*] is the vector to be solved and *<sup>I</sup>M*×<sup>1</sup> is [1, 1, . . . , 1]T.

Reference [27] has proposed the critical point of EM wave propagation. When the values of the critical points are only 0 or 1, the maximum observability can be achieved. The critical points divide the whole system into several sections. *R* detectable sections are produced; then, the propagation path of the EM wave can be obtained by the topological structure and section length of the GIS. Therefore, a constraint can be found to be

$$\mathbf{G\_{2R \times N}} \mathbf{X\_{N \times 1}} \ge \mathbf{I\_{2R \times 1}} \tag{3}$$

where *G* is the propagation path matrix of EM wave of each section and the constraint shows that at least one path of the EM wave to the node where the sensor is installed passes through both ends of the section and *<sup>I</sup>*2*R*×<sup>1</sup> = [1, 1, . . . , 1]T.

By Equations (2) and (3), the 0-1 optimization model can be obtained:

$$\begin{array}{c} \min \, \mathbf{W}\_{1 \times N}^{T} \mathbf{X}\_{N \times 1} \\ \text{s.t. } \mathbf{V}\_{M \times N} \mathbf{X}\_{N \times 1} \ge I\_{M \times 1} \\ \mathbf{G}\_{2R \times N} \mathbf{X}\_{N \times 1} \ge I\_{2R \times 1} \\ \mathbf{B}\_{2N \times N} \mathbf{X}\_{N \times 1} = \mathbf{b}\_{2N \times 1} \end{array} \tag{4}$$

where *W*<sup>T</sup> = [*w*1, *w*2,... , *wN*] is a weight vector and represents the tendency of each node to install a sensor. The elements values of *W*<sup>T</sup> are between 0 and 1. The equation constrains whether if a sensor should be installed at a certain node. If it should be installed at the node *k*, *B*(*k*, *k*) = 1 and *b*(*k*)=1, otherwise *B*(*k* + *N*, *k*) = 1. The remaining elements in *B* and *b* are 0. The model conforms to the standard form of the 0-1 programming model and it can be solved by using *brintprog* tool in MATLAB.

#### *2.3. Initial PD Location Method*

The time difference of arriving location method for GIS is based on the arriving time of the EM signals detected by UHF sensors installed at both ends of the PD section. Considering the outer radius of the cavity of GIS is much smaller than its total length, it can be assumed that the EM wave will travel along the path as the pattern in Figure 3 after PD occurs.

**Figure 3.** Location method based on time difference of arriving (TDOA).

For Figure 3, *vi* = *vj* = *c* = 0.3 m/ns is the propagation velocity of the EM wave, and the fault distance from the PD position to sensor 1 is given by

$$L\_{fi} = \frac{1}{2} [L\_{i\bar{j}} + (t\_{\bar{i}} - t\_{\bar{j}})c] \tag{5}$$

$$L\_{f\bar{j}} = L\_{i\bar{j}} - L\_{f\bar{i}} \tag{6}$$

where *Lfi* and *Lfj* are the distance from the PD position to *I* and to *j*, *ti* and *tj* are the time when sensor *I* and sensor *j* detect the signals. Equation (5) is a conventional PD location method, and the PD signals can propagate to other sensors through *I* and *j*. In view of the integrity of GIS, the propagation velocity of EM wave can be considered as remaining constant. This method lays a foundation for accurate PD locating using multiple sensors in GIS. Figure 3 can be expanded to Figure 4 below, and the EM waves travel through *I* and *j* to adjacent sensor *x* and *y*. Replacing *I* and *j* with *x* and *y* in Equations (5) and (6), the distance of extended TDOA can be expressed as

$$L\_{fx} = \frac{1}{2} [L\_{xy} + (t\_x - t\_y)c] \tag{7}$$

$$L\_{fy} = L\_{xy} - L\_{fx} \tag{8}$$

**Figure 4.** Extended TDOA location method.

Similarly, taking *x* and *j* or *i* and *y* as the ends of section, the corresponding PD position can be obtained. The extended TDOA method is the basic principle for PD locating of complex GIS. The accurate timing of optical fiber connection and the synchronization of different sampling units to picoseconds are the basis of PD locating among multiple sensors.

By the optimization Equation (4), it is assumed that *P* UHF sensors are installed at the nodes *S* = [*S*1*S*2*S*<sup>3</sup> ... *SP*] of a GIS. After PD occurs, the arriving time of the EM waves detected by the UHF sensors is *TM* = [*T*1*T*2*T*<sup>3</sup> ... *TP*] <sup>T</sup> and the minimum time in *TM* is *T*min. For a node *Sm* in *S*, all adjacent nodes are *S<sup>V</sup>* = [*SvSv2Sv3* ... *Svw*], where *w* is the number of adjacent nodes. In Figure 5, *Svr* can be any node in *SV*.

The UHF sensors at *Sm* and *Sp* form a time difference location combination when PD occurs between *i* and *j*, and the shortest path between the two nodes must pass through an adjacent section *L* of *Sm*. The other end of *L* is the node *Svr*. For *Sm* and its adjacent node matrix *SV*, *w* TDOA location combinations can be obtained. For the UHF sensor arrangement matrix *S*, *P* – 1 combinations can be obtained. Therefore, the approximate solution can be obtained by comparing the similarity between the transmission time of the EM wave from the PD position to the sensors and the arriving time of the EM wave acquired by sensors synchronously.

The Euclidean distance is used to compare the similarity between two vectors. The Euclidean distance can magnify the error caused by sensor communication and data processing. It is suitable for comparing the unknown quantity with the known quantity with possible error and detecting redundancy. The value of the Euclidean distance is inversely proportional to the similarity between the two sets of data. The more similar the two sets of data are, the smaller the value is. Therefore, the initial locating results of PD can be obtained.

**Figure 5.** The sensors adjacent to Sm.

By the shortest path algorithm, the time matrix *D<sup>f</sup>* of the EM waves propagating from the *w* × (*P* – 1) positions to UHF sensors can be found as

$$D\_f = \begin{array}{ccccc} & 1 & 2 & \cdots & w \\ & 2 & & & & \\ & & \begin{pmatrix} & \mathcal{D}\_{11} & & \mathcal{D}\_{12} & \cdots & \mathcal{D}\_{1w} \\ & \mathcal{D}\_{21} & & \mathcal{D}\_{22} & \cdots & \mathcal{D}\_{2w} \\ & \vdots & \vdots & \ddots & \vdots \\ & \mathcal{D}\_{(P-1)1} & \mathcal{D}\_{(P-1)2} & \cdots & \mathcal{D}\_{(P-1)w} \end{pmatrix} \end{array} \tag{9}$$

where the matrix element *Dij* = [*D*1 *D*2 ... *DP* ' ]. Furthermore, the distance measure matrix *E* can be built as

*E* = 1 2 ··· *w* 1 2 . . . *P* − 1 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ *E*11 *E*12 ··· *E*1*w E*21 *E*22 ··· *E*2*w* . . . . . . ... . . . *<sup>E</sup>*(*P*−1)<sup>1</sup> *<sup>E</sup>*(*P*−1)<sup>2</sup> ··· *<sup>E</sup>*(*P*−1)*<sup>w</sup>* ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ = 1 2 ··· *w* 1 2 . . . *P* − 1 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ *D*11 − *T*M *D*12−*TM* ··· *D*1*w*−*TM D*21−*TM D*22−*TM* ··· *D*2*w*−*TM* . . . . . . ... . . . *D*(*P*−1)1−*TM D*(*P*−1)2−*TM* ··· *D*(*P*−1)*w*−*TM* ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (10)

The initial PD location result is given by the minimum value in the measure matrix *E*, and the corresponding PD position can be expressed as *f* .

#### *2.4. Accurate PD Location Method for GIS*

The accuracy of PD locating can be improved by using the information of all UHF sensors installed in GIS. From the initial PD position and the inherent topological structure of GIS, the time of PD occurrence can be calculated. Then, the time of the EM wave propagating to UHF sensors is calculated by the time of PD occurrence and the EM wave velocity to obtain multiple sets of PD location results. With these results, accurate PD locating can be achieved.

The occurring time *t*<sup>0</sup> of PD can be found to be

$$T\min - t\_0 = \frac{L\_{bf'}}{c} \tag{11}$$

where *Lbf'* is the distance from the PD position to the UHF sensor corresponding to *Tmin*. By *t*0, the distance matrix of the PD signals propagating to each sensor can be calculated to be

$$\begin{array}{rcl} \mathbf{L} & = [ \ \mathbf{L}\_{1f} \ \ \ \ \mathbf{L}\_{2f} \ \ \cdots \ \ \ \ \ \ \ \ \ \ \ \ \mathbf{L}\_{Pf} \ \ ]^\mathbf{T} \\ & = \mathbf{c}(T\_M - t\_0 \boldsymbol{\eta}) \end{array} \tag{12}$$

where *<sup>η</sup>* = [1 1 ··· 1]<sup>T</sup> <sup>P</sup>×1. According to the topological structure of GIS, the PD positions can be found:

$$\mathbf{S} = \begin{bmatrix} \ S\_1 & \ S\_2 & \cdots & \ S\_P \end{bmatrix}^T \tag{13}$$

In combination with all the information of UHF sensors, the accurate PD position can be found as follows:

$$S\_{act} = \frac{1}{P} \sum\_{i=1}^{P} S\_i \tag{14}$$

The flow chart of the PD location method for GIS and the steps are as follows (shown in Figure 6):


**Figure 6.** Flow chart of the location method.

#### *2.5. Noise Reduction Method*

An improved dual-tree complex wavelet transform (DTCWT) method is proposed to reduce noise. The DTCWT method is used to decompose the noisy signal into a series of detail components which can reflect PD [28]. However, only part of DTCWT components details are closely related to PD, while the others are reflected in noise signals. Therefore, the improved DTCWT method is proposed to select sensitive components according to time-domain kurtosis and envelope spectral kurtosis. According to the principle of maximum kurtosis deconvolution, the sensitive components are deconvoluted to reduce noise by MKD deconvolution. Finally, the denoised components are reconstructed by inverse DTCWT transform to obtain the denoised signal.

The steps of the improved DTCWT method are as follows:


$$Dn = \begin{pmatrix} D1 \\ D2 \\ \vdots \\ Dn \end{pmatrix} \begin{pmatrix} d\_1(1) & d\_1(2) & \cdots & d\_1(k) \\ d\_2(1) & d\_2(2) & \cdots & d\_2(k) \\ \vdots & \vdots & & \vdots \\ d\_n(1) & d\_n(2) & \cdots & d\_n(k) \end{pmatrix} \tag{15}$$

3) Screen Signal Components. The time kurtosis of signal *Di* can be found:

$$Kur(d) = \frac{\frac{1}{n}\sum \left(d\_i - \overline{d}\right)^4}{Std(x)^4} \tag{16}$$

The envelope spectral kurtosis can be found:

$$
\omega n \upsilon(t) = \sqrt{\mathbf{x}(t)^2 + \left\{HT[\mathbf{x}(t)]\right\}^2} \tag{17}
$$

where *HT* is Hilbert transform. The higher time kurtosis and envelope spectral kurtosis indicate that there are more partial discharge impulse components in *Di*. The product of time domain kurtosis and envelope spectral kurtosis (index TE) of each component is calculated, and the component with large product is selected for reconstruction.

4) The selected DTCWT detail signal component is processed by MKD denoising, and then the inverse DTCWT transform is used to reconstruct the signal. After the inverse transform, the denoised signal of the component can be obtained.

#### **3. Experimental Results and Analysis**

#### *3.1. Experimental Platform*

The experiment is based on the partial discharge online monitoring system of China Electric Power Research Institute. The precise timing system with optical fiber connection can ensure that the clocks of each sensor can be synchronized to the same nanosecond to realize a precise location system composed of multiple channels. The test system consists of 500 kV GIS model, power supply, PD detector, phase acquisition device, built-in UHF sensor, UHF amplifier, broadband oscilloscope, discharge model, etc (shown in Figures 7 and 8).

**Figure 7.** Partial discharge (PD) online monitoring system diagram.

(**a**)

**Figure 8.** Laboratory for the PD online monitoring system. (**a**) Experimental Platform, (**b**) ultra-high frequency (UHF) sensor in GIS, (**c**) discharge model in GIS.

The whole experimental GIS platform is divided into independent test chambers by pelvic insulators. Operating hand holes are opened on the side of each chamber to facilitate the placement of defect models, such as a pin-to-plate discharge model.

#### *3.2. Noise Reduction*

The noise reduction method is validated by a measured partial discharge signal with noise shown in Figure 9.

**Figure 9.** Measured partial discharge signal.

The improved DTCWT method is used to denoise the measured partial discharge signal. The sampling number of UHF signal is 3000. The decomposition level of DTCWT is 5. The signal components of DTCWT, of which the TE values are in the top three, are reconstructed. The noise reduction result is shown in Figure 10.

**Figure 10.** Noise reduction result.

Variation trend parameter (*VTP*) and signal-to-noise ratio (*SNR*) are used to evaluate the denoising effect of UHF PD signals:

$$SNR = 20 \lg \frac{\max\_{i=1}^{n} s(i)}{\max\_{i=1}^{n} n(i)} \tag{18}$$

$$VTP = \frac{\sum\_{i=2}^{n} \left[ f(i) - f(i-1) \right]}{2\sum\_{i=2}^{n} \left[ s(i) - s(i-1) \right]} + \frac{\sum\_{i=2}^{n} \left[ f(i-1) - f(i) \right]}{2\sum\_{i=2}^{n} \left[ s(i-1) - s(i) \right]} \tag{19}$$

The closer the *VTP* value is to 1, the better the waveform fitting effect is. Comparing the DTCWT method with improved DTCWT method, the result is shown in Table 1. The improved DTCWT method has higher *SNR* and the waveform is not obviously distorted, which maintains the characteristics of the original UHF PD signal.


**Table 1.** The comparison of output variation trend parameter (*VTP*) and signal-to-noise ratio (*SNR*).

#### *3.3. Optimal Placement*

The 500 kV GIS experimental platform was used to verify the effectiveness of the proposed method. The nodes of GIS are numbered and the node number is also the sensor number. The lengths of the sections between the nodes are shown in Table 2.


**Table 2.** GIS topological parameters.

The undirected graph of the structure of this GIS shown in Figure 11 can be drawn from Table 2. The nodes where UHF sensors should be installed can be obtained by model (4), including nodes 1, 3, 4, 6, 7, 9, 11, 13, 14, 15, 16, 18, 21, 23, 25, 26, and 28. In Figure 11, the system has two tangent ring structures; therefore, non-detection zones may exist. The minimum and maximum values of the EM wave critical points for all sections are 0 and 1 through calculations [27]. Therefore, these sections are totally detectable. In addition, the constraints of model (4) ensure that the sensors can detect the EM waves effectively.

**Figure 11.** Sensor placement results.

#### *3.4. PD Location Results*

Assume a PD occurred at 20 ns in section 8-9, and the distance from the PD position to node 8 is 1000 mm. Each sensor can monitor the signal of each node and 400 MHz high-pass filtering was performed on the EM signal detected by the sensors. The amplitude and waveform of the signal barely change through the filter, so it is easy to obtain the arriving time of the EM wave, as shown in Table 3. From Table 2, *T*min = *t*7. Calculating the TDOA of sensor 7 with the other 16 sensors in sequence, the results are shown in Table 4. In Table 4, the minimum value in the Euclidean distance measure matrix corresponds to the result calculated by sensor 7 and 11. Therefore, the PD occurred in section 8-9, 954.3 mm away from node 8. The actual discharge distance is 1000 mm, the absolute error is 45.7 mm, and the relative error is 4.6%. The occurring time *t*<sup>0</sup> of PD was found to be 20.1 ns by Formula (12). Based on Formulas (13) and (14), the revised PD position can be obtained shown in Table 5. By Formula (15), the PD occurred in section 8-9, 993.5 mm away from node 8. The actual discharge distance is 1000 mm, the absolute error is 6.5 mm, and the relative error is 0.7%.



**Table 4.** Initial results of PD locating.



**Table 5.** Revised results of PD locating.

#### *3.5. Fault Tolerance Analysis*

The proposed method has been verified in the example above on PD of 500 kV GIS, but the PD locating is under the premise that there is no error in the detection of the UHF sensors. In practical operation, error data of sensors would be caused by the influence of environmental noise and communication delay.

In order to verify the fault tolerance of the algorithm, the arriving time of the EM wave of the sensor at node 11 was artificially set as 53.8 ns. The minimum arriving time *T*min remains *t*7, and the TDOA calculations of sensor 7 with the other sensors in sequence were carried out. The result shows that the PD distance calculated by the combination of senor 7 and 11 is infinity. Therefore, the information obtained by sensor 11 was judged to be erroneous data. After the erroneous data was eliminated, it can be obtained that the PD occurred on section 8-9, 988.7 mm away from node 8. The absolute error is 11.3 mm and the relative error is 1.1%.

In the same case, the existing methods for PD locating cannot eliminate the impact of erroneous data. The proposed method realizes the identification of the erroneous data through the information of multiple sensors in GIS and ensures the fault tolerance of the algorithm.

#### *3.6. Comparison of Methods*

The received signal strength indicator (RSSI) location method is based on the strength of EM wave. PD locating can be achieved by the reconstruction algorithm of RSSI based on pattern matching [29,30]. Twelve experiments were carried out to compare the proposed method and the RSSI method. Location results are shown in Table 6.


**Table 6.** Method comparison.

For Table 6, The average absolute error of the proposed method is 22.8 mm, which is obviously superior to the RSSI method, and its 100% absolute errors are less than 50 mm. The proposed method

keeps the ns-level time synchronization between the sensors and the acquisition system. Therefore, the precision and accuracy can be guaranteed.

#### **4. Conclusions**

Considering the attenuation of the EM wave when passing through the insulators, L-type and T-type structures in GIS, and the constraints of economy and observability, an optimal placement model of UHF sensors is proposed. The model can meet the requirements and has high reliability, as validated by experimental results.

Based on the optimal placement of the UHF sensors and the arriving time of EM waves, the initial results of PD locating can be calculated by the Euclidean distance and the TDOA method. Then, the accurate results can be calculated by the occurring time of PD and the information of all UHF sensors.

Compared to the existing PD location methods, the proposed method ensures the accuracy and fault tolerance of PD locating. The experimental results show the method has high accuracy and robustness.

**Author Contributions:** Conceptualization, R.L.; methodology, S.W.; validation, P.C.; formal analysis, N.P. and S.W.; data curation, Y.L.; writing—original draft preparation, S.W.; writing—review and editing, R.L.

**Funding:** This work was supported by the Fundamental Research Funds for the Central Universities (2017XKQY033). **Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
