Multi-Period linearized optimal power flow model incorporating transmission losses and thyristor controlled series compensators

ISSN 1859-1531 - THE UNIVERSITY OF DANANG, JOURNAL OF SCIENCE AND TECHNOLOGY, NO. 6(127).2018 31 MULTI-PERIOD LINEARIzED OPTIMAL POWER FLOW MODEL INCORPORATING TRANSMISSION LOSSES AND THYRISTOR CONTROLLED SERIES COMPENSATORS Pham Nang Van1, Le Thi Minh Chau1, Pham Thu Tra My2, Pham Xuan Giap2, Ha Duy Duc2, Tran Manh Tri2 1Hanoi University of Science and Technology (HUST); van.phamnang@hust.edu.vn, chau.lethiminh@hust.edu.vn 2Student at Department of Electric Power Systems, Hanoi University of Science and Technology (HUST) Abstract - This paper presents multi-period linearized optimal power flow (MPLOPF) with the consideration of transmission network losses and Thyristor Controlled Series Compensators (TCSC). The transmission losses are represented using piecewise linear approximation based on line flows. In addition, the nonlinearity due to the impedance variation of transmission line with TCSC is linearized deploying the big-M based complementary constraints. The proposed model in this paper is evaluated using PJM 5-bus test system. The impact of a variety of factors, for instance, the number of linear blocks, the location of TCSC and the ramp rate constraints on the power output and locational marginal price (LMP) is also analyzed using this proposed model. Key words - Multi-period linearized optimal power flow (MPLOPF); mixed-integer linear programming (MILP); transmission losses; Thyristor Controlled Series Compensators (TCSC); big-M 1. Introduction Electricity networks around the world are experiencing extensive change in both operation and infrastructure due to the electricity market liberalization and our increased focus on eco-friendly generation. Managing and operating power systems with considerable penetration of renewable energy sources (RES) is an enormous challenge and many approaches are applied to cope with RES integration, mainly the management of intermittency. In addition to increasing power reserves, energy storage systems (ESS) can be invested to mitigate the uncertainty of RES. The increasing application of ESS as well as problems including time-coupled formulations such as power grid planning, N-1 secure dispatch and optimal reserve allocation for outage scenarios have led to extended optimal power flow (OPF) model referred to as multi- period OPF problems (MPOPF) [1]-[2]. Typically, the MPOPF problem is approximated using the DC due to its convexity, robustness and speed in the electricity market calculation [3]. To improve the accuracy of the MPOPF model, transmission power losses have been integrated. This is significant because the losses typically account for 3% to 5% of total system load [4]. When power losses are incorporated in the MPOPF model, this model becomes nonlinear. To address the nonlinearity, reference [3] deploys the iterative algorithm based on the concept of fictitious nodal demand (FND). The disadvantage of this approach is that the MPOPF problem must be iteratively solved. Reference [5] presents another approach in which branch losses are linearized. The branch losses can be expressed as the difference between node phase angles or line flows [4]. The main drawback of this model is that it can lead to “artificial losses” without introducing binary variables [5]. Moreover, the TCSC is increasingly leveraged in power systems to improve power transfer limits, to enhance power system stability, to reduce congestion in power market operations and to decrease power losses in the grid [6]. When integrating TCSC in the MPOPF problem, this model becomes nonlinear and non-convex since the TCSC reactance becomes a variable to be found [7]. At present, there are several strong solvers like CONOPT, KNITRO for solving this nonlinear optimization problem [8]. However, directly solving nonlinear optimization problems cannot guarantee the global optimal solution. References [9]-[10] demonstrate the relaxation technique to solve the nonlinear optimization problem in power system expansion planning considering TCSC investment. Furthermore, the iterative method is used to determine optimal parameter of TCSC in reference [11]. The main contributions of the paper are as follows: - Combining different linearized techniques to convert the nonlinear MPOPF to the mixed-integer linear MPOPF. - Analysing the impact of some factors such as the number of loss linear segments, the location of TCSC as well as the ramp rate of the units on the locational marginal price (LMP) and generation output. The next sections of the article are organized as follows. In section 2, the authors present general mathematical formulation of multi-period optimal power flow (MPOPF) model incorporating losses and TCSC. The different linearization techniques are specifically presented in section 3 and 4. Section 5 demonstrates multi-period linearized optimal power flow (MPLOPF) model. The simulation results, numerical analyses of PJM 5-bus system are given in section 6. Section 7 provides some concluding remarks. 2. General mathematical formulation For normal operation conditions, the node voltage can be assumed to be flat. A multi-period optimal power flow (MPOPF) considering network constraints can be modeled for all hour t, all buses n, all generators i, and all lines (s, r) as follows: ( ) ( ) ( ), min , . , i gi gi P t T i I b G t b t P b t        (1) Subject to ( ) ( ) ( ) ( ) ( ) : , : , , 0 , g d gi dj n i i n M j j n M P t P t P t n N t T    − − =       (2) ( ) ( )  ( )max , ; , ; , ,ub lsr rs srP t P t P s r t T       (3) ( ) ( ) ( )0 , , ; , ,ubgi gi iP b t P b t i I b G t t T        (4) 32 Pham Nang Van, Le Thi Minh Chau, Pham Thu Tra My, Pham Xuan Giap, Ha Duy Duc, Tran Manh Tri ( ) ; ,lb ubgi gi giP P t P i I t T      (5) ( ) ( )1 ; ,upgi gi iP t P t R i I t T− −      (6) ( ) ( )1 ; ,dngi gi iP t P t R i I t T− −      (7) The objective function in (1) represents the total system cost in T hours (here, T = 24 h). The constraints (2) enforce the power balance at every node and every hour. The constraints (3) enforce the line flow limits at every hour. The constraints (4) and (5) are operating constraints that specify that a generator’s power output as well as power output of each energy block must be within a certain range. The other constraints included in the formulation above are the ramp-up constraints (6) and ramp-down constraints (7). If the reactance of branch xsr is taken as a variable due to TCSC installation, in the range of min max[ , ]sr srx x , it yields a new model: ( ) ( ) ( ), , min , . , sr i gi gi P x t T i I b G t b t P b t        (8) Subject to min max sr sr srx x x  (9) ( ) ( )2 7− (10) The above general model is nonlinear. Sections 3 and 4 present different linearization methods to convert this model to the linear form. 3. Linearization of the network losses In this section, the subscript t is dropped for notational simplicity. However, it could appear in every variable and constraint. Additionally, the expressions presented below apply to every transmission line; therefore, the indication ( ), ls r  will be explicitly omitted. The real power flows in the line (s, r) determined at bus s and r, respectively, are given by ( ) ( ) ( ), 1 cos sinsr s r sr s r sr s rP G B      = − − − −  (11) ( ) ( ) ( ), 1 cos sinrs s r sr s r sr s rP G B      = − − + −  (12) The real power loss in the line (s, r), ( ),losssr s rP   can be attained as follows: ( ) ( ) ( ) ( ) 2 , , ,losssr s r sr s r rs s r sr s rP P P G       = +  − (13) In the lossless DC model, the real power flow in the line (s, r) at bus s is approximately calculated as in (14): ( ) ( ) ( ) 1 ,sr s r sr s r s r sr F B X       − − = − (14) Substituting (14) in (13), the real power loss in the line (s, r) is expressed as in (15): ( ) ( ) ( ) 2 2 2 , 1 / loss sr sr s r sr sr sr sr sr sr R P G X F F R X   = = + (15) Equation (15) can be further simplified. The resistance Rsr is usually much smaller than its reactance Xsr, particularly in high voltage lines. Consequently, (15) can be further reduced to (16) ( ) 2,losssr s r sr srP R F  = (16) The first advantage of (16) compared to (13) is that power flows in lines neither built nor operative are zero. Another advantage of (16) is its possible application to model losses in HVDC lines. The quadratic losses function (16) can be expressed using piecewise linear approximation according to absolute value of the line flow variable as follows: ( ) ( ) ( ) 1 , L loss sr s r sr sr sr l P R l F l   = =  (17) To complete the piecewise linearization of the power flows and line loss, the following constraints are necessary to enforce adjacency blocks: ( ) ( )max. ; 1,..., 1sr sr srl p F l l L   = − (18) ( ) ( ) max1 . ; 2,...,sr sr srF l l p l L −  = (19) ( ) ( )1 ; 2,..., 1sr srl l l L  − = − (20) ( ) 0; 1,...,srF l l L = (21) ( )  0;1 ; 1,..., 1sr l l L  = − (22) Constraints (18) and (19) set the upper limit of the contribution of each branch flow block to the total power flow in line (s, r). This contribution is non-negative, which is expressed in (21) and limited upper by max /ubsr srp P L = , the “length” of each segment of line flow (18). A set of binary variables ( )sr l is deployed to guarantee that the linear blocks on the left will always be filled up first; therefore, this model eliminates the fictitious losses. Finally, constraints (22) state that the variables ( )sr l are binary. A linear expression of the absolute value in (17) is needed, which is obtained by means of the following substitutions: sr sr srF F F + −= + (23) sr sr srF F F + −= − (24) ( )0 1 ubsr sr srF P −  − (25) 0 ubsr sr srF P +  (26) In (24), two slack variables srF + and srF − are used to replace Fsr. Constraints (25) and (26) with binary variable θsr ensure that the right-hand side of (23) equals its left-hand side. Moreover, the slopes of the blocks of line flow ( )sr l for all transmission lines can be given by Eq. (27). ( ) ( ) max2 1sr srl l p = −  (27) It is emphasized that the number of linear segments will radically affect the accuracy of the optimal problem solution. Moreover, this linear technique is independent of the reference bus selection and thereby eliminating discrimination in the electricity market operation. Using the above expressions, the real power flow in line (s, r) computed at bus s and r can be recast as follows, respectively: ISSN 1859-1531 - THE UNIVERSITY OF DANANG, JOURNAL OF SCIENCE AND TECHNOLOGY, NO. 6(127).2018 33 ( ) ( ) ( ) ( ) 1 1 , , 2 1 2 loss sr s r sr s r sr L sr sr sr sr l P P F R l F l F      = = + = + (28) ( ) ( ) ( ) ( ) 1 1 , , 2 1 2 loss rs s r sr s r sr L sr sr sr sr l P P F R l F l F      = = − = − (29) The power withdrawn into a node n, ( ),nP t can be written as ( ) ( ) 1:( , ) 1 2l L n nk nk nk nk lk n k P R l F l F =   = +      (30) A linear substitution for the function in (3) can be found by the following equivalent constraints without increasing the number of rows ( ) ( ) 1 1 2 L ub sr sr sr sr sr l R l F l F P = +  (31) Rewriting Eq. (31), the constraints (3) are expressed as follows ( ) ( ) 1 1 1 2 L ub sr sr sr sr l R l F l P =   +      (32) 4. Linearization of a bilinear function When xsr is taken as a variable, constraint (14) also makes the MPOPF model nonlinear since this constraint is a bilinear function. To overcome the nonlinearity of this constraint, we introduce a new variable Fsr, instead of variable xsr. After obtaining the optimal solution with variable (P, F, δ), the optimal reactance can be uniquely determined according to Eq. (33) s rsr sr x F  − = (33) Therefore, the constraint (9) becomes: min maxs rsr sr sr sr x x x F  −  =  (34) It is noted that the sign of Fsr cannot be determined beforehand. Moreover, if the denominator Fsr is zero, the numerator s r − must be zero. As a result, (34) can be converted into the expression (35) depending on the sign of Fsr. min max max min 0 . 0 0 0 . sr sr sr s r sr sr sr s r sr sr sr s r sr sr if F F x F x if F if F F x F x          −   = − =    −  (35) These condition constraints can be combined by leveraging binary variables ysr and big-M based complementary constraints as follows [12]. In our model, M is taken to be / 2 due to system stability requirement [13]. ( ) ( ) min max max min1 1 sr sr sr s r sr sr sr sr sr sr s r sr sr sr My F x F x My M y F x F x M y     − +  −  +  − − +  −  + − (36) It is important to stress that linear technique using the above binary variable is exact while the linearized technique in Section 3 is approximately presented. 5. Multi-period linearized optimal power flow (MPLOPF) model with losses and TCSC The MPLOPF model with losses and TCSC has the following form: ( ) ( ) ( ), , min , . , i gi gi P F t T i I l G t b t P b t        (37) Subject to ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) : , : , 1 :( , ) 1 1 1 , , 2 ; , , , g d l gi dj i i n M j j n M L nk nk nk nk l L L k n k nk nk l l P t P t R l F l t F l t n t F l t F l t    + − =  + − = = − =    +           + −             (38) ( ) ( ) ( ) 1 1 1 , , 2 L ub sr sr sr sr sr l R l F l t F l t P + − =    + +       (39) ( ) ( ) ( )max, . , F , ; 1,..., 1sr sr sr srl t p F l t l t l L + −  + = − (40) ( ) ( ) ( ) max, , 1, . ; 2,...,sr sr sr srF l t F l t l t p l L + −+  −  = (41) ( ) ( ), 1, ; 2,..., 1;sr srl t l t l L  − = − (42) ( ) ( ) ( )  , 0; , 0; , 0;1sr srF l t F l t l t + −  = (43) ( ) ( ) ( ) 1 0 , ; , , L ub l sr sr l F l t t P s r t T+ =       (44) ( ) ( ) ( ) 1 0 , 1 ; , , L ub l sr sr l F l t t P s r t T− =   −       (45) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) min max max min 1 1 sr sr sr s r s r sr sr sr sr sr sr s r s r sr sr sr My t F t x t t t t F t x My t M y t F t x t t t t F t x M y t         − +  −   −  +  − − +  −     −  + −    (46) ( ) ( )4 7− (47) Regarding the computational complexity of the model, the number of continuous variable is 24. .GEN GENiN N ( )24. 1 2.24. .BUS LINN N L+ − + and the number of binary variables is ( )24. . 1 2.24.LIN LINN L N− + . After the MPLOPF problem is solved, the marginal cost at the node i in hour t can be determined by the following expression [3]: . .i E i E l i l l LMP LMP LF LMP SF −= − +  (48) 6. Results and discussions In this section, the multi-period linearized optimal power flow model is performed on the modified PJM 5-bus system [3]. The MPLOPF problem is solved by CPLEX 12.7 [15] under MATLAB environment. 34 Pham Nang Van, Le Thi Minh Chau, Pham Thu Tra My, Pham Xuan Giap, Ha Duy Duc, Tran Manh Tri 6.1. System data The test system is shown in Figure 1. The total peak demand in this system is 1080 MW and the total load is equally distributed among buses B, C and D. The daily load curve is depicted in Figure 2. Two small size generators on bus A have the capability to quickly start up. The ramp rate for the other generators is 50% of the rated power output [14]. E D A B C Limit=240 MW Brighton Park City Load Center Solitude Sundance 110MW $14 600MW $10 200MW $35 520MW $30 100MW $15 Figure 1. PJM 5-bus system and generation parameters Figure 2. Daily load curve for PJM system 6.2. Impact from the number of linear blocks Table 1. The effects of number of linear blocks Linear blocks Objective ($) Total losses (MW) Time (s) 2 3844.43 316.69 1.71 4 3824.04 244.83 2.97 6 3822.96 238.56 5.28 8 3820.70 230.41 8.42 10 3820.55 229.49 12.35 11 3820.51 229.49 14.61 The number of linear blocks can significantly affect the solution time as well as the model accuracy listed in Table 1. The key idea in this paper is to find the number of linear blocks which give the best balance between the model accuracy and the solution time. In this case, 10 is an appropriate number in terms of objective value, total losses and calcultaion time. 6.3. Impact from losses Table 2 compares the results of power output at 10 AM using the proposed model. These results are also compared with those of POWERWORLD software using the ACOPF model [16]. When comparing to POWERWORLD software, the calculated results using the proposed model considering losses are more accurate and less different than that of the model neglecting losses. Table 2. Generating output results at 10 AM Bus Lossless (MW) Losses (MW) POWERWORLD (MW) A1 110 110 110 A2 100 100 100 C 19.95 30.1 27.83 D 195.05 194.8 197.2 E 600 600 600 Figure 3. LMP at bus B at different hours without losses and with losses The results of LMP calculations at node B for 24 hours using the proposed model with and without losses are given in Figure 3. This figure illustrates that the effect of power losses on LMP is very little. This result is consistent because the power losses account for about 1% of the total load for this PJM 5-bus system, therefore the marginal generating units as well as congested lines are the same in both cases. 6.4. Impact from TCSC location It is assumed that power losses are not considered and the ramp rate of the generating units (not including units at node A) are taken as 25% of the maximum power output. Also, the compensation level of TCSC varies from 30% to 70%. Figure 4 depicts the power output of generator at node C for 24 hours for different locations of TCSC. During the period from 1 AM to 3 AM, the power output of the unit at node C nearly remains when the location of TCSC varies. In addition, the power output of this unit is highest in 24 hours when TCSC is located in line A-B. Figure 4. The dependence of Generating output of Unit at bus C on TCSC location 6.5. Impact from ramp rate constraints Figure 5 shows the power output of generator located at node C when changing the ramp rate of generators and it is assumed that TCSC is not applied to the power grid. From the 5 AM to 24 PM, the power output of this unit is the same for ramp rates of 25%, 35% and 50%. At the same time, the output of this unit is the highest for ramp rate 100% of the maximum power. 900 950 1000 1050 1100 0 5 10 15 20 25 Lo ad ( M W ) Hour 20 25 30 35 0 5 10 15 20 25 L M P ( $ /M W h ) Hour Losses Lossless 0 200 400 0 5 10 15 20 25 G en er a ti o n ( M W )) Hour Line A-B Line B-C ISSN 1859-1531 - THE UNIVERSITY OF DANANG, JOURNAL OF SCIENCE AND TECHNOLOGY, NO. 6(127).2018 35 Figure 6 depicts the effect of TCSC placement on the power output with different ramp rate scenarios at 10 AM. We see that the power output of generator at node C does not change as the ramp rate of the units changes in case of placing TCSC on line A-B. However, when TCSC is not installed, the ramp rate of units has a significant effect on the unit's output, increasing from 30,097 MW for the ramp rate of 50% to 223,37 MW for the ramp rate of 100%. Thus, using TCSC also reduces the impact of the ramp rate on the power output. Figure 5. The dependence of generating output of Unit at bus C on Ramp rate without TCSC Figure 6. The dependence of power output of Unit at bus C on Ramp rate with TCSC in line A-B at 10 AM 7. Conclusion This paper presents multi-period linearized optimal power flow (MPLOPF) model based mixed-integer linear programming (MILP). This MPLOPF integrates line losses and Thyristor Controlled Series Compensator (TCSC). The different linearization techniques, such as piecewise linear approximation and big-M based complementary constraints are deployed to convert multi-period nonlinear OPF problem to multi-period linearized OPF model. The calculated results using the proposed model are compared to those of the commercial POWERWORLD software and this proves the validation of the proposed model. Additionally, the influences of the number of linear blocks, line losses, location of TCSC and ramp rate are analyzed. The results reveal that these factors can importantly impact on LMP, generating output of units as well as revenue of participants in electricity markets. NOMENCLATURE The main mathematical symbols used throughout this paper are classified below. Constants: ( )sr l Slope of the lth segment of the linearized power flow in line (s, r) ( ),gi b t Offered price of the bth linear block of the energy bid by the ith generating unit in hour t srB Imaginary part of the admittance of line (s, r) srG Real part of the admittance of line (s, r) srR Resistance of the line (s, r) srX Reactance of the line (s, r) ( )djP t Power consumed by the jth load in hour t L Number of the blocks of the loss linearization ub srP Transmission limit of line (s, r) ub giP Upper bound on the power output of the ith producer lb giP Lower bound on the power output of the ith producer up iR Ramp-up limit of the ith unit dn iR Ramp-down limit of the ith unit min srx Lower bound of the reactance of the line with TCSC max srx Upper bound of the reactance of the line with TCSC BUSN Number of nodes GENN Number of generators LINN Number of transmission lines GEN iN Number of energy blocks of unit i Variables: ( ),giP b t Power output corresponding to the bth block of the ith unit in hour t ( ),nP t Power withdrawal at bus n in hour t ( ),srP t Power flow in line (s, r) at node s in hour t ( ),rsP t Power flow in line (s, r) at node r in hour t ( )s t Voltage angle at node s in hour t ( )srF t Power flow in line (s, r) in hour t without losses ( ),losssrP t Power losses in line (s, r) in hour t ( )sr l Binary variable relating to the line flow linearization ( )sry t Binary variable corresponding the big-M based complementary constraints ( )srx t The reactance of the line with TCSC in hour t iLF Loss factor at bus i l iSF − Sensitivity of branch power flow l with respect to injected power i l Shadow price of transmission constraint on line l Sets: I Set of indices of the generating units ( )iG t Set of blocks energy bid offered by the ith unit in hour t N Set of indices of the network nodes l Set of transmission lines ACKNOWLEDGMENT This research is funded by the Hanoi University of Science and Technology (HUST) under project number T2017-PC-093. REFERENCES [1] D. Kourounis, A. Fuchs, and O. Schenk, “Towards the next generation of multiperiod optimal power flow solvers”, IEEE Trans. 0 200 400 1 3 5 7 9 11 13 15 17 19 21 23 G en er a ti o n ( M W ) Hour Ramp rate 25% Ramp rate 35% Ramp rate 50% 0 50 100 150 200 250 25% 35% 50% 100% G en er at io n ( M W ) Ramp-rate TCSC in line A-B No TCSC 36 Pham Nang Van, Le Thi Minh Chau, Pham Thu Tra My, Pham Xuan Giap, Ha Duy Duc, Tran Manh Tri Power Syst., vol. 8950, pp. 1–10, 2018. [2] P. N. Van, N. D. Huy, N. Van Duong, and N. T. Huu, “A tool for unit commitment schedule in day-ahead pool based electricity markets”, J. Sci. Technol. Univ. Danang, vol. 6, pp. 21–25, 2016. [3] F. Li, S. Member, R. Bo, and S. Member, “DCOPF-Based LMP Simulation : Algorithm, comparison with ACOPF and sensitivity”, IEEE Trans. Power Syst., vol. 22, no. 4, pp. 1475–1485, 2007. [4] D. Z. Fitiwi, L. Olmos, M. Rivier, F. de Cuadra, and I. J. Pérez- Arriaga, “Finding a representative network losses model for large- scale transmission expansion planning with renewable energy sources”, Energy, vol. 101, pp. 343–358, 2016. [5] J. M. Arroyo and A. J. Conejo, “Network-constrained Multiperiod auction for a pool-based electricity market”, IEEE Trans. Power Syst., vol. 17, no. 4, pp. 1225–1231, 2002. [6] P. N. Van, N. D. Hung, and N. D. Huy, “The impact of TCSC on transmission costs in wholesale power markets considering bilateral transactions and active power reserves”, J. Sci. Technol. Univ. Danang, vol. 12, pp. 24–28, 2016. [7] G. Y. Yang, G. Hovland, R. Majumder, and Z. Y. Dong, “TCSC allocation based on line flow based equations via mixed-integer programming”, IEEE Trans. Power Syst., vol. 22, no. 4, pp. 2262–2269, 2007. [8] Alireza Soroudi, Power System Optimization Modeling in GAMS. Springer, 2017. [9] O. Ziaee, O. Alizadeh Mousavi, and F. Choobineh, “Co- optimization of transmission expansion planning and TCSC placement considering the correlation between wind and demand scenarios”, IEEE Trans. Power Syst., vol. 8950, no. c, pp. 1–1, 2017. [10] M. Farivar and S. H. Low, “Branch Flow Model: Relaxations and Convexification (Parts I, II)”, pp. 1–11, 2012. [11] P. N. Van and L. M. Khanh, “The optimal location and compensation level of Thyristor Controlled Series Compensator (TCSC) in Wholesale Electricity Markets considering Active Power Reserves”, J. Sci. Technol. Tech. Univ. Vietnam, 2017. [12] T. Ding, R. Bo, W. Gu, and H. Sun, “Big-M Based MIQP Method for Economic Dispatch With Disjoint Prohibited Zones”, IEEE Trans. Power Syst., vol. 29, no. 2, pp. 976–977, 2014. [13] T. Ding, R. Bo, F. Li, and H. Sun, “Optimal Power Flow with the Consideration of Flexible Transmission Line Impedance”, IEEE Trans. Power Syst., vol. 31, no. 2, pp. 1655–1656, 2016. [14] Y. Wei, H. Cui, X. Fang, and F. Li, “Strategic scheduling of energy storage for load serving entities in locational marginal pricing market”, IET Gener. Transm. Distrib., vol. 10, no. 5, 2016. [15] IBM, “IBM ILOG CPLEX Optimization Studio Community Edition”. [16] https://www.powerworld.com/ (The Board of Editors received the paper on 18/4/2018, its review was completed on 04/5/2018)