HOW TO EFFICIENTLY INCORPORATE FACTS DEVICES IN OPTIMAL ACTIVE POWER FLOW MODEL

. This paper presents for the ﬁrst time how to easily incorporate facts devices in an optimal active power ﬂow model such that an eﬃcient interior-point method may be applied. The optimal active power ﬂow model is based on a network ﬂow approach instead of the traditional nodal formulation that allows the use of an eﬃciently predictor-corrector interior point method speed up by sparsity exploitation. The mathematical equivalence between the network ﬂow and the nodal models is addressed, as well as the computational advantages of the former considering the solution by interior point methods. The adequacy of the network ﬂow model for representing facts devices is presented and illustrated on a small 5-bus system. The model was implemented using Matlab and its performance was evaluated with the 3,397-bus and 4,075-branch Brazilian power system which show the robustness and eﬃciency of the formulation proposed. The numerical results also indicate an eﬃcient tool for optimal active power ﬂow that is suitable for incorporating facts devices.


1.
Introduction. The concept of using solid state power electronic converters for power flow control at the transmission level has been know as facts (Flexible AC Transmission Systems). Facts devices have a large potential ability to make power systems operate in a more flexible, secure and economic way. The possibility of operating power systems at the minimal cost while satisfying specified transmission and security constraints is one of main current issues in stretching transmission capacity by the use of facts devices. The conventional optimal power flow model, however, must undergo some changes such as the inclusion of new control variables representing facts devices and the corresponding adaptation of the load flow solution.
Since the 70s, some optimal active power flow problems have been formulated using network flow models as an alternative to the classical approach based on the nodal formulation [4,16,5,8,11,9,7,19,10,1,3,14,6,12]. Several applications have been suggested such as reliability studies [16], load supply capability [5], economic dispatch [8,12], contingency analysis [7], hydrothermal scheduling [3], available transfer capability [18], and transmission network planning [11,9,19]. One advantage of the network flow approach is that power flows are explicitly represented in the model. This feature allows transmission capacities being handled directly as limits imposed on the transmission variables, the consideration of power flows dependent functions, such as losses, as part of the optimization criteria [9,10,12], and suitable representation of facts devices [6,20].
The main purpose of this paper is to present an efficient active optimal power flow model that easily incorporates and handles facts devices. This is important since the objective is solve large scale power system problems in an affordable time and also present solutions that considers facts equipments. A case study with the Brazilian power system, composed of 4075 branches and 3397 buses, is selected for performance evaluation of the optimal active power flow model and the computational savings due to sparsity exploitation are highlighted. The approach is based on a network flow model and is solved by an interior point method (IPM) [12]. It is shown that the IPM computational effort does not increase with the presence of facts devices, as occurs with the traditional nodal formulation [13,17], because the nodal formulation, when include facts equipments, do not allow the use of matrix sparsity exploitation in an interior point method.
The modeling of facts devices, such as variable phase shifter (VPS), thyristorcontrolled series compensator (TCSC) and thyristor-controlled phase shifter (TCPS), are addressed in the paper. A 5-bus system illustrates the incorporation of these devices in the network flow model.
The outline of the paper is as follows. Notation is presented in Section 2. The network flow model for the optimal active power flow problem, and its equivalence with the nodal formulation, are presented in Section 3. Extensions of the network flow model to incorporate facts devices are described in Section 4. In Section 5 some computational results for the Brazilian power system are presented. Finally, Section 6 summarizes the conclusions.

Notation.
m: number of branches. n: number of bars. l: number of loops. g: number of generators.
A: network incidence matrix (n × m). L: network loop matrix (l × m). X: reactance diagonal matrix (m × m). R: resistance diagonal matrix (m × m). B: susceptance matrix (n × n). E: matrix (n × g) formed by columns with all elements equal zeros except the element of a line corresponding to a generating bar. p: active power generation vector (g × 1). d: active power demand vector (n × 1). f : active power flow vector (m × 1). ϕ: phase shifter angle vector (m × 1). θ: angular vector (n × 1). f min : lower bound for active power flow f . f max : upper bound for active power flow f . p min : lower bound for active power generation p. p max : upper bound for active power generation p. α,β: weight constants. φ 1 : function associated with power flow vector. φ 2 : function associated with generation vector. * : symbol for fixed or target value.
3. Optimal active power flow modeling. In an active power flow model, Kirchhoff's Current Law (KCL) states power flow conservation on the buses, and can be expressed by Kirchhoff's Voltage Law (KVL), by its turn, states that branch power flows are proportional to angular differences and susceptances, and can be expressed by The classical nodal formulation of the active power flow problem merges both laws into a single linear equation by substituting Equation (2) into Equation (1) where B = A(X −1 A t θ) is the susceptance matrix. Therefore, the active power flow model based on nodal formulation reduces the linear system that represents Kirchhoff's laws to a dimension of the number of busses by eliminating branch power flow variables and expressing Kirchhoff's laws only in terms of power injections and bus angle variables. The nodal formulation brings considerable reduction on the dimension of the Kirchhoff's law linear system, from (n + m) constraints and (2n + m) variables in Equations (1) and (2) to (n) constraints and (n + m) variables in Equation (3). For this reason, it has been widely adopted specially in applications where branch power flows or its limits does not need to be known.
For applications where branch power flows have to be taken into consideration, additional linear constraints should be considered for branch power flow variables being calculated from bus angle variables through Equation (2). In other words, KVL is considered twice since it is implicitly present in Equation (3) and explicitly present in the additional linear constraint of Equation (2).
In order to overcome the difficult of handling power flow transmission variables and constraints through the nodal formulation of the active power flow model, an usual approach is to apply a relaxation procedure by which only the transmission limits which are violated are progressively enforced through dual simplex iterations. Of course this approach introduces a considerable additional computational effort since it requires an interactive solution procedure as well the monitoring of transmission power flow violations [15].
The indirect handling of transmission variables in the nodal formulation of the active power flow problem also constitutes a difficulty for representing facts devices. By this reason, one suggestion is to return the original Kirchhoff's Laws given by Equations (1)-(2) [17].
In another approach, followed in this paper, the active power flow problem is formulated by a network flow approach where the Kirchhoff's laws are expressed by the following set of linear equations since LA t θ means that angular differences along all basic loops must be zero and Equation (5) express this in terms of the active power flow variable f [1]. The explicit use of f variable, as used in Eq. (5), in the formulation of KVL is an advantage when is necessary to consider power flow limits in optimal active power flow model. Another advantage of the Eq. (5) over the Eq. (2) is the reduction on the number of equations necessary to represent KVL since the first depends on the number of basic loops l = n − m + 1 and the second depends on the number of branches (m). Therefore, the network flow approach given by Equations (4)-(5) is more compact than the original Equations (1)-(2) since it considers KVL only in the branches that belong to loops. Efficient procedures to find the loop matrix L from the incidence matrix A are discussed in [12,2], being the approach described in [12] adopted in this paper. It is important to note that, although it could be laborious to identify the set of loop equations in large scale systems, this task must just be done once.
Once the network flow is proved to be a more compact model to represent Kirchhoff's laws, another important step is the construction of an active optimal power flow model by defining the representation of two subjects: objective function and bounds in power flow and generation variables.
The optimal active power flow problem can then be formulated as the following bounded network flow model with additional linear constraints The objective function (7) is the composition of two different criteria, the first one depending on power flows, φ 1 (f ), and the second one depending on power generations, φ 2 (p). Both criteria are represented by quadratic and separable functions, and are combined through scalar weights α and β in a bi-objective optimization framework.
φ 1 (f ) can be a quadratic separable function expressed by: where M 1 , m 2 , m 3 are given diagonal matrix, vector and scalar parameters, respectively. By setting M 1 = R, m 2 = 0 and m 3 = 0, function φ 1 (f ) can represent an approximation of the transmission power losses in the network. φ 2 (p) can be a quadratic separable function expressed by: where N 1 , n 2 and n 3 are given diagonal matrix, vector and scalar parameters, respectively. By setting appropriated values for N 1 , n 2 and n 3 , function φ 2 (p) can represent quadratic generation costs. A quite useful objective function φ 2 (p) is the quadratic deviation from a desirable generation dispatch. Such a dispatch can arise from a pool auction in an electricity market or from a dispatch model which does not take into consideration transmission network constraints. In such case φ 2 (p) can be represented by where W is a diagonal matrix with the component w i as the penalty term associated with deviation from the desired generation p * i . The Equation (14) could be derived from Eq (13) by setting N 1 = W , n 2 = −p * W and n 3 = 1 2 p * W p * . Note that model (7)-(11) corresponds to a bounded nonlinear minimum cost network flow problem with additional linear constraints, whose special structure allows efficient solution by IPM techniques. The recently predictor-corrector IPM proposed by [12] efficiently solves the optimal active power flow problem modeled by the network flow approach. A key factor for computational performance of the solver is to observe that major computational burden is due to solving the following linear system  ( where K = A LX , D f is a proper diagonal matrix related with power flows and r is a proper vector. The addition of any canonical column e j from E to K in order to get a square nonsingular matrix was proposed in [12]. Using this feature with the Sherman-Morrison-Woodbury formula an efficient solution of (15) was achieved. A 5-bus system illustrates the network flow model. Figure 1 presents the singleline diagram and Figure 2 presents the corresponding equivalent graph with two basic loops.
For the graph of Figure 2, and setting all decision variables in the left-hand-side, constraints (8) and (9) can be rewritten as: In the traditional nodal formulation of the optimal active power flow problem, bus angles are assigned as decision variables instead of power flows, resulting the following model where V = X −1 A t and B = A(X −1 )A t . Note that Equation (17) represents both Kirchhoff's laws through a single constraint since it results from substituting Equation (2) into Equation (1), that is Note also that for fixed susceptances, Equation (17) imposes limits on power flows through limits on the corresponding angular differences. Since the directly handling of such constraints would carry out a great computational burden, the usual nodal formulation handles transmission limits through an interactive relaxation procedure by which current violated transmission limits are gradually imposed by dual simplex iterations [15].
Both the network flow model (7)-(11) and the nodal model (17)-(17) can be expressed through a mathematical programming problem of minimizing a quadratic separable objective function subject to linear constraints and bounds.
M in z t Hz + g t z + c st : where the meaning of variables and parameters are given on Table 1.
As can be seen, the network flow model is less compact than the nodal model, which seems to be an disadvantage. Both models have power generations as variables (g), being the difference that the network flow model have also branch power flows (m) as variables whereas the nodal model have bus angles (n). The number of branches (m) is greater than the number of buses (n) by the exact amount of basic loops in the network (l = m − n). Thus, the network flow model has the same number of additional variables (m − n) and constraints (l) with respect to the nodal model. For the Brazilian electrical network with 4075 branches and 3397 buses, the increase in system dimension is around 20 per cent for the general model proposed by Eq. (19). Nevertheless, system structure on the network flow model can be better exploited by the interior-point method (IPM) than the nodal one, and this is an aspect that can greatly overcome the dimension disadvantage created by the general model (Eq. (19)) if adequately exploited.
For the network model (Eqs. (7)-(11)) the IPM can be directly applied, because all matricial constraints must be equalities. But to prepare the nodal model (Eqs. (17)-(17)) for the IPM development is necessary to transform the Eq. (17) in a equality constraint, resulting in the formulation (20) described in [13]: This new nodal formulation has (m + n) constraints and (m + n + g) variables, while network model has less constraints (m) and less variables (g + m), and the two models have g degrees of freedom.
An improvement in the nodal model can be made if the second equality constraint is applied in the first constraint, giving a more simplified nodal model as described by [17]: This nodal model has (m + n) constraints and (m + n + g) variables, but does not need the calculation of B = AX −1 A t matrix for the KCL constraint.
The formulation (21) can be transformed in the network model (Eqs. (7)- (11)) if the property LA t = 0 is applied in the equation V θ = f . Also note that this will eliminate the variable θ from the network flow model resulting in a more compact model with less variables (g + m) and constraints (m).
There are some limitations for the development of an efficient IPM that complains the two nodal models exposed.
First, the nodal model showed in (20) deals with matrix B = AX −1 A t by removing the line (and column) corresponding to the slack bar, and is less advantage on sparsity because B is more dense in comparison with K, which preserves the network structure of the problem expressed by matrices A and L.
For the nodal model showed in (21), the branch limits could be handled by two approaches : to explicit f as a variable in the model or to consider the introduction of matrix inequality constraint. Both have problems. In other words: • The model as given by (21) has three groups of variables and this will destroy matrix sparsity structure proposed by [12], because it will produce the following K matrix: • The second equality constraint and the bounds on the f variable could be merged in the inequality f min ≤ V θ ≤ f max . Since the efficient IPM proposed by [12] complains only equality constraints, the inequality will produce slack variables which will destroy the matrix sparsity structure and consequent exploitation.

4.1.
Modeling variable phase shifters. For simplicity, it is first considered a fixed phase shifter device. In this case, constraints (8) and (9) should be modified to Considering the 5-bus system of Figure 1, and placing a fixed phase shifter device in branch 2 − 5, results the following linear system As can be seen, a simple change in the right-hand-side of Equation (9) can cope with the presence of fixed phase shifters in the network.
To include the phase shifter as a variable, which occurs when an specific power flow value should be assigned to a transmission line, the linear system turns to be as follows In order to assign a fixed power flow value f * 25 at branch 2 − 5 of the 5-bus system, a variable phase shifter is placed in that branch, resulting the following linear system Note that a simple change of the columns corresponding to the envolved variables (power flows and phase shifter angles) can cope with the consideration of variable phase shifter devices in the network.

4.2.
Modeling TCSC facts. The thyristor-controlled series compensator (TCSC) consists of a capacitor bank series shunted by a thyristor-controlled reactor to provide a smoothly variable capacitive reactance series. Its consideration in the network will modify the linear system as follows Again, note that the modifications of the linear system to consider TCSC devices are quite simple and will not lead to increase in computational burden.

4.3.
Modeling TCPS facts. The controllable parameter of thyristor-controlled phase shifter (TCPS) is the voltage shift angle Ψ that produces a active power injection in the bars of the branch with this equipment as follows where Ψ N are the variable voltage shift angle for the branches with TCPS device.
To inclusion of TCPS in branch 2 − 5 will produce the follow matricial result:  (22) can be built from A and LX matrices through droppings and permutations of columns, accordingly to the specific device modeling and without any additional computational burden. This property remains correct even two or more of devices are considered in the same network, because if one parameter in a branch is fixed then another parameter must be variable and viceversa. In other words there is a kind of duality between the parameters.
This means that the inclusion of power flow control devices will not change the number of variables and constraints associated with Kirchhoff's laws. In particular, observe that duality between the parameters is responsible to keep KVL always linear.
Moreover, the network flow model with facts devices will result in a sparsity structure that still permits the speedup granted by the use of Eq. (15).
This advantage is not shared by the Nodal model (17)- (17), because the consideration of facts devices introduce additional constraints [13], or turn the it nonlinear [17]. 5. Numerical results. In order to evaluate the efficiency of the network flow model for optimal active power flow problems, the Brazilian Power System (BPS) with 4075 branches and 3397 bars was selected.
The first case study presents the active power flow problem with the objective of minimizing the quadratic deviation from a given operation point, considering light, medium and heavy loads. In a first case, both the network flow and the traditional nodal models were considered without sparsity exploitation. In a second case, the network flow model with sparsity exploitation was considered. Table 2 shows the CPU time and the number of iterations of each case study. Note that both models have similar CPU time per iteration and number of iterations, which leads to similar total CPU time. However, the network flow approach with sparsity exploitation speed up the CPU time in almost 70 times. These savings are not possible with the nodal model due to the dense structure of the corresponding linear system. A second case study was performed considering different setups for a TCSC device with the objective of minimizing transmission losses in the BEP system under heavy load. In this case study different fixed values for the power flow in the TCSC branch were considered, and the minimal loss dispatch and the corresponding TCSC reactance values were determined by the network flow model.
Although the BPS have more than 150 generation bars, the impact of different TCSC setups carried out changes in generation of only six plants, as show in Fig.  3 and Fig. 4.
The most significant change happens in the generations of Tucuruí, Lajeado and Serra da Mesa plants when the flow in TCSC FACT is negative (North is sending energy although the rain sazonality indicates a contrary flow orientation). Figure (4) gives a geographical insight of the results for different fixed values of power flow in the TCSC device. This study is carried when the plants of the north region which Tucuruí belongs have a low volume and its necessary to import energy from the southwest system which Lajeado and Serra da Mesa belongs. This TCSC equipment have been installed for this situation and its natural flow orientation is to be a positive one.
Expanding this study and obtaining the resulting transmission loss and the TCSC reactance values for more fixed values of flow in the TCSC FACT for the high load level scenario in BIN produces the Fig. 4.
The Figure 5 indicates that the optimal flow value for the TCSC in heavy load level in BPS is around the base value f * (f * = 84.5 MW, is the point named). Another important phenomena occurs when the flow in the TCSC is negative, and then the reactance of TCSC will be put in its maximum value (−0.0054 pu), denoting that the flow in the TCSC will be controlled only by generation changes.
This explains why in the case with f = −300 happens a significant change in generation of Tucuruí (generation increase) and S.Mesa (generation reduction) in comparison with the case with no flow between north and southwest regions (f = 0).  6. Conclusion. This paper has presented for the first time how to efficiently incorporate facts devices in an optimal active power flow model solved by an interiorpoint method. One advantage of the proposed methodology is that it conjugates an interior point method with network flow modeling which allows the use sparsity exploiting and reduce the computational time 70 times when compared with a nodal approach (which cannot use sparsity exploitation in an interior-point method approach). The presented model is also robust and a quadratic separable objective functions of active power generation and flows also can be considered, thus allowing the minimization of active transmission power loss, generation cost or quadratic deviation from a desired dispatch, or a combination of both. The model was tested on three load scenarios for the 3,397-bus and 4,075-branch Brazilian Power System which show the robustness and efficiency of the formulation proposed. Another advantage is that the network flow allows the incorporation of FACTS equipments without any impact on the computational time to solve the problem. The results showed that the technique is flexible, robust and efficient and the adequacy of the interior-point method to incorporate FACTS devices