Analysis of Optimal Steady-State Operation of Power Systems with Embedded FACTS Devices: A Matlab-Based Flexible Approach

This book chapter presents a flexible approach to incorporate mathematical models of FACTS devices into the Power Flow (PF) and the Optimal Power Flow (OPF) analysis tools, as well as into the standard OPF Market-Clearing (OPF-MC) procedure. The proposed approach uses the Matlab Optimization Toolbox because it allows to easily: (a) implement a given optimization model, (b) include different objective functions using distinct equality and inequality constraints and (c) modify and reuse an optimization model that has been previously implemented. The conventional OPF model is the main core of the proposed approach, which is easily implemented and adapted to include the mathematical models of FACTS devices. The resulting implementation of the OPF model featuring FACTS devices can be easily modified and adjusted to obtain the implementation of both the PF and the OPF-MC models which includes such devices. It should be mentioned that with the flexible approach proposed here, the complexity as well as the implementation time of optimized models featuring embedded FACTS devices is significantly reduced, since it is not necessary to define the expressions associated with the hessian matrix and the gradient vector. The flexibility and reliability of the proposed approach are demonstrated by means of several study cases using test as well as real power systems.


Introduction
The efficient operation and economics of electric power systems have always occupied an important place in the electric power industry [1]. Here, it is worth mentioning that highly desirable operative conditions cannot be successfully achieved by means of using the old electromechanically controlled devices, such as the mechanical phase shifter. These old devices execute control actions through mechanical actuators, and thus provide a poor high-speed control. Some of the consequences of this lack of speed and reliable control are associated with stability problems. Also, the power flow deviates through secondary transmission lines that are not the intended ones. In addition to the above mentioned issues, there are instability issues when the transmission resources are closely used to their thermal and economic limits. Not to mention, the high or low voltages resulting from such a poor high-speed control, in addition to many other issues [2]. Those poor operating conditions have given place for the need to develop faster controllers and a more efficient power system management. Fortunately, the fast development of power electronics based on new and powerful semiconductor devices has given rise to novel technologies such as FACTS controllers that greatly improve the operating conditions of power systems [3].
From a steady-state operation perspective, the control capabilities of FACTS devices allow to adjust the active and reactive power flows at their output terminals, as well as to achieve the local control of reactive power and voltages at the connected nodes [4]. In order to investigate the effectiveness of these control capabilities for alleviating a large variety of problems associated with the steadystate performance of power systems, the existing tools for the steady-state condition assessment have been upgraded to incorporate steady-state models of FACTS controllers. The Power Flow, the Optimal Power Flow and the OPF Market-Clearing procedures are the common tools used to assess the steady-state operation and market-clearing conditions of power systems. Accordingly, the upgrade of those three tools has received great attention. For instance, references [4-12, 4, 13-19, 20-24] are representative samples of the several publications where models of FACTS controllers are incorporated into the Power Flow, the Optimal Power Flow and the market-clearing procedures, respectively.
Certainly, the evolution of FACTS controllers is going to be progressive as function of time [2]. With the further development of technologies based on power electronics, the economic viability of these controllers will improve in such a way that more and more FACTS devices are expected to be designed and used in future applications [13]. In this sense, it is expected that modern power systems are going to be electronically controlled rather than mechanically controlled [2]. Consequently, to study the usefulness of new FACTS devices for improving the steady-state performance of power systems as well as to speed up research efforts, it would be needed the upgrade of the existing steady-state analysis tools in order to incorporate the corresponding mathematical models of such devices.
Bearing in mind the aforementioned prospective on FACTS controllers, in this book chapter is proposed a Matlab-based flexible approach that incorporates mathematical models of FACTS devices into: the power flow and the optimal power flow analysis tools; as well as into the standard OPF market-clearing procedure [25]. The flexibility of the proposed approach is achieved by taking advantage of both: (a) the tractability of the nonlinear continuous optimization theory to consider different objective functions that can be properly constrained and (b) the simplicity of the Matlab Optimization Toolbox to implement and solve an optimization model, as well as its ability to modify and reuse an optimization model that has already been implemented [26]. In order to show the flexibility of the proposed approach, first consider that the OPF model is the main core of this proposal, then the strategy to code the objective function, the equality and inequality constraints of the conventional OPF model are given in detail. Second, the implemented conventional OPF model is enlarged to include the steady-state models of FACTS controllers. Lastly, the resulting OPF implementation is modified and reused, to straightforwardly obtain the computational implementation of both the Power Flow and the standard OPF market-clearing models that use embedded FACTS controllers.
In this way, the proposed book chapter can reduce the computational implementation load of the FACTS model devices into one of the three previously mentioned procedures. Therefore, this proposal facilitates carrying out research in order to assess the effectiveness of FACTS devices for improving the steady-state performance of power systems. Arguably, the use of high-level programming languages and optimization toolboxes such as the ones provided by Matlab could reduce the computational efficiency of any approach that is implemented using such tools. However, for research purposes, the computational efficiency of the proposed approach is quite reasonable, since the solution of the distinct optimization models related to medium scale power systems is commonly obtained in a time lapse that varies from a few seconds to a maximum of 10 minutes. In addition, when the main objective is to investigate the effectiveness of FACTS devices for improving the steady-state performance of power systems, the computational efficiency requirement can be further relaxed.
This book chapter is organized as follows. Since the conventional OPF model is the core of this approach, its formulation is given in Section 2, as well as the steadystate models of the traditional power system components. These models are then considered in Section 3 to assemble the conventional OPF model. The Section 4 shows the incorporation of the FACTS models into the conventional OPF model. This last model can be straightforwardly extended in order to obtain the computational implementation of the PF and, the standard OPF market-clearing models in Sections 5 and 6, respectively. The applicability of the proposed approach for the steady-state operation of power systems using embedded FACTS devices is assessed by means of numerical examples in Section 7. Last but not least, the conclusions of this chapter book are stated in Section 8.

The optimal power flow general-formulation
The optimal power flow formulation is a non-linear optimization problem which consists in minimizing an objective function subject to a finite set of equality and inequality constraints. Mathematically, the OPF model is given by

Minimize fy ðÞ
(1) where f (y) is the objective function to be minimized, h(y) and g(y) are a finite set of functions corresponding to the equality and inequality constraints, respectively. Whilst y, represents the full set of variables of the power system with a lower limit y min and an upper limit y max . It must be mentioned that the set y contains of both the state variables as well as the control variables.
The OPF formulation is given in polar coordinates taking in consideration the model of each power system component. The most common components of a power system are: generators, loads, fixed shunt compensation elements, transmission lines and transformers; whose steady-state models are given in [4]. These models are adopted in the flexible approach proposed in the present book chapter.

Flexible approach to implement the OPF model
This section describes the flexible approach to implement the conventional OPF model as a computer program, which is developed based on structured programming along with the use of the fmincon optimization function of Matlab [26]. This proposal allows to obtain a computational implementation of the OPF model with the following remarkable features: (a) It considerably simplifies the computational complexity of the implementation of a conventional OPF, (b) it is very flexible when considering diverse physical laws and operating conditions, (c) it allows for the integration of FACTS models and (d) it has the tractability of easily being adapted to study electric power systems by means of different steady-state analyses. For example, it can be straightforwardly adapted to solve the conventional power flow problem and to carry out the analysis of electric power markets. In this way, steady-state analyses of distinct power systems featuring different electric components can be carried out by doing very few changes to the composition of the set of steady-state variables as well as the set of constraints and the objective function of the proposed implementation of the OPF model. The aforementioned features are further exploited in Sections 5, 6 and 7 in order to take into consideration: FACTS devices, the conventional power flow problem and the standard OPF market-clearing procedure, respectively.

The fmincon function of the Matlab optimization toolbox
The fmincon function of the Matlab Optimization Toolbox is used to solve the general OPF model which is a constrained nonlinear optimization or programming problem. That Matlab function attempts to find out the constrained minimum of a scalar objective function of several variables, starting with an initial estimation of the constrained minimal point. Here it is worth mentioning that the fmincon function is designed to solve problems where the objective and constraint functions are both continuous as well as their first derivatives. In order to obtain the optimal solution of the optimization problem at hand, the fmincon users can choose one of the following optimization options: the trust-region-reflective method, the active-set method, the interior-point method and the sequential quadratic programming approach [26]. Figure 1 shows a schematic representation of the main input and output arguments of this function. It must be pointed out that, depending on to the requirements of the specific optimization problem, some empty input arguments must be provided. The next command line syntax of the fmincon function, including input and output arguments, is shown below [26]. Observe that the first line of code corresponds to the input arguments, whilst the second one corresponds to the output arguments.
x = fmincon(fun,x0,A,b,Aeq,beq,lb.,ub,nonlcon,options) [x,fval,exitflag,output,lambda,grad,hessian] = fmincon(___) In Figure 1 and the syntax, the input argument fun is the handle of the M-function file containing the objective function to be optimized. X 0 represents the initial conditions of the state vector X of the optimization model, which is given by (5). For the purpose of this work, this state vector (X 0 ) stores the initial conditions of the steadystate variables of the power system. By the way, the lower and upper bounds of X 0 are assigned to the arguments lb and up, respectively. The size of lb and up must be equal to the size of the vector X 0 . Next, the nonlcon input argument is the handle of the Mfunction file containing the nonlinear inequality as well as the equality constraints. In the same venue, the argument options is a structure where the user can define optional parameters for the optimization process, such as: the optimization algorithm, the convergence tolerance, the maximum amount of iterations, the tolerance on the constraint violation and the termination tolerance on the function value, among many other parameters. It is worth mentioning that the options structure has the parameter LargeScale that allows for a choice on which algorithm to use. If that parameter is set to on, then the gradient of the objective function as well as the constraints must be provided by the user. It must be highlighted, however, that in this work, the LargeScale parameter is set to off with the aim of simplifying and making more flexible the proposed computational implementation. Last but not least, the input arguments A, b, Aeq and beq are associated with the coefficients of the linear inequality and equality constraints. In this proposal, however, the OPF model only considers nonlinear equality and inequality constraints and, therefore, the above mentioned input arguments are considered as empty arrays.
Next, the output arguments are briefly described to provide more clarity on the approach proposed here. First, X* stores the optimal value of the optimization model variables as they are assessed by the fmincon function. Second, the fval output returns the value of the objective function given by the optimal solution X*. Third, the output argument is a structure that contains information about the results of the optimization. Forth, lambda is an argument storing the value of the Lagrange multipliers at X*. Fifth, the arguments hessian and grad store the value of the Hessian matrix and the gradient vector of the objective function at X*, respectively. Last but not least, exitflag is an integer identifying the termination condition of the optimization algorithm; that is, if (a) exitflag > 0: the algorithm reached the optimal solution X*, if (b) exitflag = 0: either the maximum amount of function evaluations or iterations specified by the user was exceeded, thus the optimal solution has not been reached and if (c) exitflag < 0: the algorithm did not reach the optimal solution, therefore, there is an unfeasible solution.
A detailed way to define the set of steady-state variables, the objective function and the constraints sets of the OPF model is given below. Also, the explicit OPF formulation is given in polar coordinates in accordance to the model of each power system component, as it is described in the following subsections.

Steady-state variables
Within the proposed implementation, the set of the steady-state variables of the OPF model is treated as a vector containing such variables in the order shown in (5). That vector is called the vector of steady-state variables, denoted by X. Note that T denotes the transpose of the vector X.
In (5), the number of nodes and generators of the power system are represented by nb and ng, respectively. θ i and V i represent the angle (grades) and magnitude (pu) of the phasor voltage in each node i of the power system, and P gi is the power of the generator i. Note that the size of the steady-state vector is (2nb + ng, 1). It should be mentioned that in the implementation, the vector X also corresponds to the vector of the initial conditions of the state variables X 0 , which is also called by the fmincon function during the optimization process. That function returns the optimal value of the steady-state variables, which is stored in the vector X, that is, the result is rewritten in X. In this way, when the optimization process is finished, X is then the vector of the optimal solution represented by X*.
The vector X 0 , which represents the initial conditions of the state variables of the power system in the OPF model, is defined as follows: the magnitude and angle of the voltage in each node of the power system are initialized at 0 grades and 1.0 pu, respectively, whilst P g0 is the initial condition of active power at the generation nodes. It is worth mentioning that P g0 is set to the values of the active power obtained from the analysis of the Optimal Generation Dispatch (OGD), which in turn, neglects the losses in the transmission elements. Last but not least, observe that the initial conditions of the active power generation help to reduce the number of iterations associated with the OPF optimization process. The syntax of the vector X containing the state variables is the following: %Phase angle of nodal voltage in radians X(nb+1:2*nb,1)=V_mag; %Magnitude of nodal voltage in pu X(2*nb+1:2*nb+ng,1)=Pg; %Active power generation

Objective function
The objective function considered in the OPF model corresponds to the total generation cost of active power, which is a nonlinear function given by where the cost curve coefficients of each generator i are represented by ai, bi and ci. ng is the number of generation nodes including the slack node. P gi is the individual generation power level. Observe that the objective function in (6) is coded into an M-function file with the name objfun_OPF. This last file is also called by the fmincon function through the input argument fun, using the syntax @(X) objfun_OPF(X). The general Matlab syntax for the objective function given by Eq. (6) is the following: function f = objfun_OPF(___) f=0; for k=1:ng f=f+ai(k,1)+bi(k,1)*Pgi(k,1)+ci(k,1)*Pgi(k,1)^2; end It is worth mentioning that in the above M-function file is possible to: include, modify, define or write the objective function as it is desired. Therefore, the objective function (6) can be modified or even changed by a new function. This characteristic provides great flexibility to the proposal, including the possibility for handling several objective functions at the same time.

Equality constraints
The equality constraints of the OPF analysis are given by the balance energy equations for the power system. These restrictions must be unconditionally met because they establish the active and reactive power balance under steady-state conditions for all the nodes of the power system. If at least one constraint is not satisfied, then the solution of the OPF problem is not feasible. That set of equality constraints is given by hy ðÞ¼ where i and k are the nodes in which the energy balance is satisfied. P i inj,j and Q k inj,j represent the active and reactive injected-power at node i and k, respectively. Note that the injection point is the j-th network element. Also note that the activepower and reactive-power (denoted by P gj and Q gj , respectively) are provided by the controllable complex power sources located at the j-th generation node (j=i, k). Last but not least, the complex power consumption is due to the active-power load P lj and the reactive-power load Q lk .
In the proposed flexible implementation, the set of equality constraints is written in an M-function file with the name constraints_OPF. This set of constrains is given by (8 In the above expression, sn is the subscript associated with the slack node, whilst X sn stores the value of phase of the slack node in the corresponding component of the state vector X. Note that θ sn is the value of the phase at which the angle of the slack node is fixed. It must be pointed out that in (7) and (8), the value of the reactive power Q gi must always be zero for the i-th generation node. This is because the level of reactive power generation depends on the state variables of the power system. Henceforth, Q gi lacks a scheduled value in the optimization model OPF. Furthermore, the constraint of reactive power balance is satisfied only at nongeneration nodes. Then such constraint should be handled as a functional inequality constraint for generation nodes in order to establish the reactive power balance in such nodes, as it is explained below. Note that the way in which the vector of equality constraints (Ce) is written, provides the flexibility to add more constraints of this type. This in turn allows extending the OPF model to other similar models for the steady-state operation assessment of power systems. The general command line syntax for the vector of equality constraints, Ce, is the following: It is important to note that in the above syntax, Pg_i, Pl_i and Pinj_i are implemented as M-function files. The same is valid for the code line that corresponds to the reactive power.

Inequality constraints
The physical and operational limits are mathematically modeled as inequality constraints. This is because these limits help to constrain the range of the practical steady-state operation of the power-system components. For a better handling of the inequality constraints, these are classified into two types: inequality constraints on variables and functional inequality constraints. The limits on the magnitude of the substation-voltage and the generation of active-power are represented by the inequality constraints given by (9).
In (9), the superscripts min and max represent the lower and upper bounds of the respective variable. As mentioned before, the lower bound is assigned to the input argument lb, and the upper bound is assigned to the input argument up; both arguments are called by fmincon function during the optimization process. The assignment of such limits to the respective previously mentioned input-arguments is shown in (10). The superscript T, on lb and up indicates the transposition of the respective vector.
Clearly, the size of the vectors lb and up is equal to the size of the vector of the steady-state variables, X. As mentioned before, storing the inequality constraintvectors into an M-file enables the convenient modification of such vectors by adding or removing constraints that are associated with other models of steadystate operation of power systems. This is shown in subsequent sections.
The level of reactive power generation is not a scheduled value because it is a function of the system variables. Therefore, the limits of the power (Q gi ) are modeled in the OPF problem as a functional inequality constraint as follows: The reactive power generation level Q gi in (12) is given by the following function: In accordance with (13), the balance of reactive power at the i-th generation node is achieved only when the level of that power is inside its own bounds. If the level of the reactive power generation of a given generator is out of its upper or lower limit, then the corresponding inequality constraint (9) must be added to the active set in (12), becoming an equality constraint that helps in resetting the reactive power Q gi at the violated limit Q g,v .
It should be noted that in the OPF model, the above expression avoids the violation of the limits of the reactive power generation. Also it allows to maintain the balance of reactive power at the generation nodes. Next observe that the functional inequality constraints are written in the same file of the equality constraints, as it is shown below. %-Qg+Qlow<=0 ci(ngen+1:2*ngen,1) = Qgi -Qg_max; % Qg-Qsup<=0 As mentioned above, the input argument nonlcon is the handle of the Mfunction file containing both the nonlinear equality and the functional inequality constraints. Therefore, the two sets of constraints are written in an M-function file whose name is constraints_OPF. The fmincon function calls this M-function file through the input argument nonlcon, using the syntax @(X)Constraints_OPF(X).
In addition to the M-functions: objfun_OPF and constraints OPF; the fmincon function also calls the vectors: lb and, ub; which contain the lower and upper limits of the state variables. The general syntax for the above is the following: The way in which the above M-files are written in Matlab facilitates the modification of the constraints, the objective function and the vector of state variables, providing flexibility to the proposed implementation. Furthermore, the development of this proposed methodology is based on structured programming, which increases its flexibility by allowing the inclusion of other electrical devices and/or different steady-state models for operational assessment of distinct power systems with minor code modifications.

Structured programming of the OPF model
The proposed flexible approach is developed using structured programming because in this way it is possible to develop computer programs that can be easily understood. Also this kind of programming is very useful when it is necessary to modify or extend an existing program. This is the case of the proposal presented here, where the OPF model is extended to consider several optimization applications of power systems with embedded FACTS controllers.
In that figure, it can be observed that the structure of the proposed implementation contains a main program. This one starts reading the M-file of the power system data, which includes a convergence tolerance, the MVA-base and the data of the most common components of a power system in accordance to Section 3 of this book chapter. After reading the M-file, in the main program, the power system data are converted to pu in order to work in the same units all entire system components, but also to avoid optimization scaling problems. Next, the program calls the respective functions for calculating the conductance and the susceptance of the transmission lines, the transformers and the shunt compensators. Observe that these quantities are used for calculating the injections of active and reactive power. This calculation involves the equality constraints associated with the power balance at each node of the power system. After that, the OGD M-function file is executed to obtain the solution of the Optimal Generation Dispatch. It must be pointed out that in the OGD M-function file, the fmincon function calls the objective function and, the OGD constraints through the argument nonlcon using a similar syntax to the one men-

Including FACTS models into the proposed flexible approach
Models of FACTS devices are easily integrated into the proposed flexible implementation by considering the following issues: the state variables and their limits; the injection of active and reactive power at the node where the controller is connected and the control equations or control modes of each FACTS device. Below, it is described in detail how the FACTS devices are integrated into the flexible approach.

State variables of the FACTS devices
The models of FACTS devices add one or more state variables to the OPF formulation; therefore, it is necessary to consider such variables in the proposed implementation. The state variables of FACTS devices are included in the set of steadystate variables (5), denoted by the state vector X. It is convenient to write these variables in the next position from P gi ,thatis,intheposition(2nb+ng+1,1)ofthe vector X, in accordance to the structure of the proposed approach. When the model of a FACTS device is incorporated, it should be noted that the size of the vector X is increased from (2nb + ng, 1) to (2nb + ng + nvar, 1), where nvar is the number of state variables of each FACTS device. Similarly, the initial conditions of the state variables of FACTS devices are included in the vector X 0 in a position following the initial condition of P gi . The operating limits of FACTS devices are handled as inequality constraints on the respective sate variables. Therefore, the lower and upper bounds of the state variables for these devices are written in the position (2nb + ng + 1, 1) of the vector lb and up, respectively. Note that, when considering the model of any FACTS device, the size of these last vectors is extended to equal the size of the state vector X. The syntax of the vector X which contains both the power system state variables and the state variables added by FACTS devices is the following: %Set of state variables added by i-th FACTS device VP(2*nb+ngen+1:2*nb+ngen+nvar,1)=FACTS_i;

Power flow injections and the control equations of FACTS devices
The power balance in all nodes of a power system must be unconditionally satisfied. Therefore, the incorporation of the model of a FACTS device into the optimization model implies to consider the injection of active and reactive power provided by such device. Thus, the set of equality constraints (7) is extended to include the power flow injections provided by FACTS devices as follows: where the terms P i inj,FACTS and Q i inj,FACTS represent the active and reactive power flow injections provided by FACTS devices. In the structure of the proposed flexible approach, shown in Figure 2, the power-flow injections are included in an M-function file whose name is P_injected. This last file is called by the M-function file Constraints_OPF, according to (8). Note that the size as well as the syntax of the vector, Ce, are not modified by the inclusion of the active and reactive powerflow injections corresponding to the FACTS devices.
In addition to the extended power-balance equality constraints set, given by (16), the inclusion of FACTS devices into the optimization model implies to include their respective control equations as elements of the following set of additional equality constraints given by where F is the parameter to be controlled by the i-th FACTS controller with the specified target value F spec . Observe that such control parameter can either be the active/reactive power flow through the transmission line or the magnitude of the voltage at a specific node of the power system in accordance to the type of the FACTS device connected to the network at such node. It is noteworthy that when a device is connected to the network, its corresponding constraint is activated in the set given by (17) in order to enable the respective control action of such device during the process of solving the optimization model. Otherwise, that is, when no FACTS devices are involved, the equality constraints of the control modes remain inactive.
Thus, in the proposed flexible implementation, the vector Ce of equality constraints is extended to include the additional set of equality constraints. Mathematically, (8) is augmented to contain the equation given by (18) Adding control equations increases the size of the vector Ce, given by (8), to (2nb + 1 + nF,1), where nF equals the number of control modes of FACTS devices in the optimization model. It is noteworthy that when the model of the FACTS device has operation equations such as in the case of the VSC-HVDC system; these equations must be considered as equality constraints associated with the corresponding device operation. Also note that such constraints are included in the optimization model given by (17), in the same way as the control constraints were included in that model. Last but not least, the vector Ce must be extended to include the above equality constraints. The control modes of the FACTS devices are implemented in M-function files in order to facilitate their handling and to provide flexibility to the proposed implementation. The general syntax for including the additional set of equality constraints associated with the control equations corresponding to the FACTS devices is the following: ce(2*nb+2:2*nb+1+nF,1)=Control_FACTS(___); %nF control modes The previously mentioned procedure can be applied to integrate easily the model of any FACTS device into the proposed flexible approach. The proposed implementation has been structured to connect and disconnect any FACTS device belonging to the power system. This is achieved by writing a 0 or 1 in the data file of the corresponding FACTS controller, where 0 is for disconnecting and 1 is for connecting.

Solution of the conventional power flow problem
The conventional Power Flow analysis is aimed to determine the steady-state operation of an electric power system. This consists on computing the magnitude and angle of each voltage at all the nodes of the system. This is carried out in addition to determining the active and reactive power flow in the transmission elements of the power system. In this analysis, four quantities are associated with each node of the system which are the magnitude of the voltage, the angle of the phase and the active and reactive power-flow. Furthermore, the system nodes can be classified as: slack nodes, regulated or PV nodes and load or PQ nodes. The complex node voltage given in phase and magnitude can be specified at any slack node. Note that in the regulated nodes, the magnitude of the voltage and active power are known. Whilst the value of the active and reactive power is defined in any load node.
The power flow analysis is based on the power balance equations at each node of the power system which are given by the set of equality constraints in (7). The proposed approach in this chapter book allows obtaining the solution to conventional Power Flow problems by performing small changes in the OPF formulation. This can be achieved by converting the OPF model into a conventional PF model. In the OPF analysis is not important to identify the node types, as it is in the case of a PF analysis. However, both analyzes are based on the power balance equations for all nodes of the power system. Thus, it is possible to include the PV nodes in the OPF model by setting the magnitude of the voltage as well as the active power to a given value, only at the nodes where a generator is connected. Similarly, the slack node can be included by setting the magnitude of the voltage and the angle of the phase to a given value. In this way, the solution of the conventional PF analysis is obtained from the OPF model. The aforementioned procedure is accomplished by including an additional set of equality constraints in the optimization model. This allows for the setting of the magnitude of the voltage as well as the active power to certain values at the generation nodes. It also allows for the setting of the magnitude of the voltage and the angle of the phase to the stablished values at the slack nodes. This additional set of equality constraints that aims to obtain the solution of the PF problem is given by (19).
In this way, the conventional Power Flow analysis is modeled as a problem of a constrained nonlinear optimization which is formulated using the OPF model. Therefore, the general model of the conventional power flow analysis considered in this work is given as follows: where the h OPF-PF is the set of equality constraints given by (19). This allows the inclusion of the node types in order to convert the OPF model into a conventional PF model. Therefore in accordance to the conventional PF model given by (20)-(24), the proposed flexible approach solution to the Power Flow problem can be obtained by considering the full of the power flow model given by (6), (7), (9), (12), and (19). It is important to note that by including the above equality constraints, the resulting OPF implementation is modified and reused to straightforwardly obtain the computational implementation of the Power Flow analysis that considers embedded FACTS controllers. In this case, it must be considered the extended optimization model with FACTS devices given by (6), (9), (12), (16), (17) and (19).

Minimize fy ðÞ
In order to maintain the structure of the proposed implementation, the last set of equality constraints is written in the vector Ce of the M-function file named Constraints_OPF as follows: where the dependent variables, such as the nodal voltage, are represented by y. Whilst the control variables, for example, the supply bids P S and power demand P D , are denoted by p. The functions f, g and h are defined below.

Objective function
The objective function of the OPF market-clearing model represents the social welfare associated with the power system. That function is given by the difference between: the sum of the accepted demand bids P Di times their corresponding bid prices C Di (in $/MWhr) and the sum of the accepted production bids P Sj times their corresponding bids prices C Sj (in $/MWhr). Mathematically, this objective function is given by In (32), the first and second term is also identified as f D and f S , respectively; nl represents the total demand bids, and ng is the total production bids.

Equality constraints
The set of equality constraints (33) represents the standard power flow equations gy , p ðÞ ¼ g θ, V, k G , P S , P D ðÞ ¼ 0 where the set y=(θ, V, k G ) and p=(P S , P D ). The variables θ and V are the complex-voltage phases and magnitudes at nodes of the power system, respectively. Whilst k G is a scalar variable used to account for the system losses by means of either a unique or a distributed slack node. It is important to note that y represents the set of dependent variables to be optimized. The generator and load powers are defined as follows: where P G0 and P L0 stand for the generator and load parameters that are not part of the market trading, that is, these parameters represent the negotiation of energy between particulars in a specific big industry and/or Generation Companies (GENCOS). Whilst PS and PD are the amounts of energy that can be traded in the electricity market by GENCOS and Energy Services Company (ESCOS). Finally, the reactive loads are assumed to have constant power factors which are computed by using the following expression: PLoad=X(2*nb+ng+1:2*nb+ng+nl,1); f=0; fd=0; for k=1:PDi fd=fd+di(k,1)+ei(k,1)*PLoad(k,1)+fi(k,1)*PLoad(k,1)^2; %CDi(PDi) end fs=0; for k=1:ng fs=fs+ai(k,1)+bi(k,1)*Pgi(k,1)+ci(k,1)*Pgi(k,1)^2; %CSj(PSGj) end f=fs-fd;

Numerical examples
Two numerical examples are tested (with the six-node power system and, an equivalent model of the Mexican Interconnected Power System -MIS.), to numerically illustrate the effectiveness and easy implementation of the proposed flexible approach to carry out the analysis of the OPF market-clearing, as well as the OPF and the conventional PF procedures. The numerical examples presented in this section were carried out in Matlab-R2018a, using a personal computer with a processor Intel Core i5-3210M CPU running at 2.50 GHz with a 6 GB of RAM.
In that way, the solution for the standard OPF market-clearing is obtained, in addition to solve both the OPF and the PF models. For both power systems, the numerical examples are designed as follows. First, the solution of standard OPF market-clearing model is obtained. Next from that solution, the value of the sum of the inelastic and the elastic loads is assigned to the corresponding load nodes in order to solve the OPF model. In other words, the Optimal Power Flow analysis is carried out by considering the value of the load obtained by OPF market-clearing study as an active load. Finally, the conventional PF analysis is executed by setting both the magnitude of the voltage and the active power generated at the PV nodes, to their corresponding values found by the OPF analysis. Similarly, within this last algorithm step, the voltage phases and magnitudes in the slack node are set to their corresponding values found by the OPF analysis. Here, it is worth mentioning that the value of the active load obtained in the market-clearing analysis is considered as a load. Also observe that the objective functions are the same in the three models mentioned above.

Testing the six nodes power system
The test case of the six node power system reported in [25] is reproduced in this section in order to take the corresponding results as a reference to determine the solution of the OPF and conventional PF models. The test power system consists of three GENCOS and ESCOS, respectively, which are interconnected through 11 transmission lines, as shown in Figure 3. The complete simulation data set for this system is given in [23,25]. For the test cases belonging to this test power system, the convergence tolerance is set to 1 Â 10 À9 .
The steady-state solutions achieved by the proposal to: the standard OPF market-clearing, the OPF and the PF models are presented in Table 1. Here, it is worth mentioning that such solutions correspond to the case of not using FACTS devices.
As expected, Table 1 shows that the solution of the three procedures is practically the same. This is because the numerical examples were designed to obtain those results. The active and reactive loads computed by the OPF market-clearing procedure for ESCO 1, 2 and 3 are: 115 MW-60 MVar, 110 MW-70 MVar and 99.767 MW-54.418 MVar, respectively. The CPU time for the case of the OPF market-clearing model is 2.3 s. Whilst for the case of the OPF and the PF models, such a time measure yields 2.1 and 1.5 s, respectively. Subsequently, the six nodes power system is modified to simultaneously include the TCSC and SVC controllers. The TCSC is commissioned to maintain the active power flow within the range of 5.62-7.0 MW across the transmission line connecting nodes 4 and 5. The first terminal of the TCSC controller is connected to node 4, whilst the second terminal is connected to a new substation called 6a. Note that the SVC is placed at node 5 to keep the magnitude of the voltage of this node at 1.02 pu. For the numerical  examples, the state variables representing the thyristor's firing angle of the TCSC and the SVC devices. These angles are both initialized at 155°with a lower and an upper limit of 90 and 180°, respectively. The fixed values of the capacitive and inductive reactances of the controllers are X C = 0.9375% and X L = 0.1625%, respectively. These values are referred to a voltage base of 400 kV and a power base of 100 MVA. Table 2 shows the corresponding steady-state solutions given by the three analyses for the case of a steady state operation using embedded FACTS devices. Clearly, the value of the objective function is increased, due to the TCSC sets the active power flow through a transmission line at higher value than in the base case. This control action causes a redistribution of power flow in the system, increasing the active power losses and the accepted demand bids. Regarding the control performed by the SVC, there is a decrease in the generation of reactive power as it is shown in Table 2. This is due to the reactive power losses are keep to a minimum limit in the power system by the SVC controller. This controller supplies reactive power in order to control the magnitude of the voltage at the substations. Therefore, such controller serves as a source of reactive power into the power grid. Although the integration of these FACTS devices into the power system increases the power generation cost and decreases the value of the social welfare net, it is very important to install them in the power system. This is because such devices can improve the steady-state operation of the power system, which by the way can be evaluated using the proposal presented here.
When the FACTS devices are embedded in the power system, the CPU time required for obtaining the solution of the OPF market-clearing model is 2.3 s. Whilst for the OPF and the PF models, the steady-state solution is determined in 2.2 and 1.6 s, respectively.

The Mexican interconnected power system (MIS)
The system considered in this numerical example is a reduced model of the Mexican Interconnected System. The MIS is represented by an equivalent model composed of: 190 nodes, 46 generators, 90 loads, and 265 transmission lines operating at voltage levels ranging from 115 to 400 kV [27].
The limits for the magnitude of the voltage for all nodes are set to 0.94 ≤ V i ≤ 1.07 pu. Whilst the active and reactive power limits for all generators are set to: 200 ≤ P gi ≤ 1000 MW and À 250 ≤ Q gi ≤ 350 MVar, respectively. The cost functions of all generators are considered linear with a value of the linear coefficient between 0.0019 and 0.0020 $/MWhr. The MIS consists of two areas identified as Area A and Area B, interconnected through a relatively weak double-circuit tie line, as it is shown in Figure 4.  DC voltage is fixed at 3.0 pu. Whilst the DC link resistance is set to 0.0034 pu. The summary results of the OPF market-clearing, the Optimal Power Flow and, the conventional Power Flow, all of them including the PtP VSC-HVDC system embedded within the MIS system, are presented in Table 3.
Observe that instead of the VSC-HVDC, the Phase-Shifting Transformer (PST) is used to maintain the active power flow through the inter-area link at the same level in MW as in the above case. Also note that the primary and secondary winding impedances of the PST exhibit zero resistance. Then the primary and secondary inductive reactances are set to 0.0 and 0.05 pu, respectively. Next, the complex tap ratios are T v = Uv = 1.0∠0°. Note that the control of the active power flow is carried out by means of a primary phase angle control. For implementing such control, the limits of the phase-shifter angle are set to AE10°. Last but not least, the summary results of the three steady-state cases previously considered in this book chapter are shown in Table 3, where a PST controller has been included.
The above table shows that when the FACTS controllers are integrated into the power system, the value of the objective function increases in the case of the OPF and the PF studies, but it decreases in the OPF market-clearing case. This is caused by the variations of the active power losses in the transmission elements, which in turn are due to the control of the active power flow performed by the VSC-HVDC and PST devices through the inter-area link. In the cases of the market-clearing studies using: (i) the PST and (ii) the VSC-HVDC; 20 and, 22 generators hit their limits of P g , respectively. Whereas for the OPF and the PF studies using: (i) the PST and (ii) the VSC-HVDC; 13 and, 10 generators hit their limits of P g , respectively.
Regarding the magnitudes of the nodal-voltages, 17 of them at powersubstations hit their limits when the analysis of the market-clearing is performed using the PST and VSC-HVDC. Whereas in the OPF and the PF analysis featuring both FACTS devices, a total of four magnitudes of the substation-voltages hit their limits. It should be mentioned that the variables that violated their limits were either fixed at their upper or lower limits in a proper manner by the flexible approach proposed here. The CPU times required to obtain the solution of the market VSC-HVDC, the market PST, the OPF VSC-HVDC, the OPF PST, the PF VSC-HVDC and the PF PST are: 1503, 978, 66, 64, 34 and 33 s, respectively. From the results of the numerical examples presented in this section, it can be inferred that the proposed implementation reliably carries out the steady-state operation assessment of power systems using embedded FACTS devices. The value of the upper row corresponds to f S and the lower row is f D . Table 3.
Summary results of the steady-state studies with the VSC-HVDC and PST devices.

Conclusions
A Matlab-based flexible approach to carry out the analysis of the steady-state operation of power systems using embedded FACTS controllers has been presented in this book chapter. The flexibility and reliability of this proposal has been demonstrated by means of several numerical examples. The referred numerical examples show that by using the proposed approach, it is possible to simulate the optimal operation of power systems in a simple way, because it avoids the implementation of complex programming algorithms typical of the majority of the optimization models that corresponds to power systems. For example, this is the case for the complexity related to the program code lines related to the hessian matrix as well as the gradient vector.
In this regard, it is also important to mention that the proposed flexible approach is easy to implement as a computer program in order to carry out the analysis of the optimal steady-state operation conditions corresponding to power systems featuring FACTS devices. This is because the proposed flexible approach greatly reduces the complexity as well as the implementation time of such optimization models.
Besides that, the flexibility of the proposed approach provides the possibility to modify and/or consider different objective functions as well as equality and inequality constraints. Therefore, in addition to solving the optimization models considered in this work, the steady-state model of any other FACTS device could also have been used. Not to mention that a large variety of steady-state optimization applications for power systems could also been readily implemented by using the flexible approach proposed here.
Last but not least, the computational efficiency of the proposal is quite reasonable at least for research and academic purposes, since most of the solutions to the optimization models related to medium scale power systems can be usually obtained within a few seconds to about 10-15 minutes. Furthermore, it is possible to infer that this proposed approach may lay a solid basis for researchers, teachers and students who seek to develop their own computational tool, whether it is for research, training or educational purposes.