PV Curves for Steady-State Security Assessment with MATLAB

Most of the problem solutions oriented to the analysis of power systems require the implementation of sophisticated algorithms which need a considerable amount of calculations that must be carried out with a digital computer. Advances in software and hardware engineering have led to the development of specialized computing tools in the area of electrical power systems which allows its efficient analysis. Most of the computational programs, if not all of them, are developed under proprietary code, in other words, the users does not have access to the source code, which limits its usage scope. These programs are considered as black boxes that users only need to feed the required input data to obtain the results without knowing anything about the details of the inner program structure. In the academic or research areas this kind of programs does not fulfill all needs that are required hence it is common the usage of programming tools oriented to the scientific computing. These tools facilitate the development of solution algorithms for any engineering problem, by taking into account the mathematical formulations which define the solution of the proposed problem. Besides, it is also common that most of these programs are known as script or interpreted languages, such as MATLAB, Python and Perl. They all have the common feature of being high level programming languages that usually make use of available efficient libraries in a straightforward way. MATLAB is considered as a programming language that has become a good option for many researchers in different science and engineering areas because of it can allow the creation, manipulation and operation of sparse or full matrices; it also allows to the user the programming of any mathematical algorithm by means of an ordered sequence of commands (code) written into an ascii file known as script files. These files are portable, i.e. they can be executed in most of software versions in any processor under the operating systems Linux or windows. The main objective of this chapter consists in presenting an efficient alternative of developing a script program in the MATLAB environment; the program can generate characteristic curves power vs. voltage (PV curves) of each node in a power system. The curves are used to analyze and evaluate the stability voltage limits in steady state, and they are calculated by employing an algorithm known as continuous load flows, which are a variation of the Newton-Raphson formulation for load flows but it avoids any possibility of singularity during the solution process under a scenario of continuous load variation. To illustrate the application of this analysis tool, the 14-node IEEE test system is used to


Introduction
Most of the problem solutions oriented to the analysis of power systems require the implementation of sophisticated algorithms which need a considerable amount of calculations that must be carried out with a digital computer.Advances in software and hardware engineering have led to the development of specialized computing tools in the area of electrical power systems which allows its efficient analysis.Most of the computational programs, if not all of them, are developed under proprietary code, in other words, the users does not have access to the source code, which limits its usage scope.These programs are considered as black boxes that users only need to feed the required input data to obtain the results without knowing anything about the details of the inner program structure.In the academic or research areas this kind of programs does not fulfill all needs that are required hence it is common the usage of programming tools oriented to the scientific computing.These tools facilitate the development of solution algorithms for any engineering problem, by taking into account the mathematical formulations which define the solution of the proposed problem.Besides, it is also common that most of these programs are known as script or interpreted languages, such as MATLAB, Python and Perl.They all have the common feature of being high level programming languages that usually make use of available efficient libraries in a straightforward way.MATLAB is considered as a programming language that has become a good option for many researchers in different science and engineering areas because of it can allow the creation, manipulation and operation of sparse or full matrices; it also allows to the user the programming of any mathematical algorithm by means of an ordered sequence of commands (code) written into an ascii file known as script files.These files are portable, i.e. they can be executed in most of software versions in any processor under the operating systems Linux or windows.The main objective of this chapter consists in presenting an efficient alternative of developing a script program in the MATLAB environment; the program can generate characteristic curves power vs. voltage (PV curves) of each node in a power system.The curves are used to analyze and evaluate the stability voltage limits in steady state, and they are calculated by employing an algorithm known as continuous load flows, which are a variation of the Newton-Raphson formulation for load flows but it avoids any possibility of singularity during the solution process under a scenario of continuous load variation.To illustrate the application of this analysis tool, the 14-node IEEE test system is used to integrating solution algorithms for economic dispatch, calculation and plotting of PV curves, testing of methods with distributed slacks, static var compensator (SVC) models, transmission lines, generators, etc.It is also an important tool in power system research, making the PTL more complete for the analysis or studies oriented to the operation and control of electric power systems.

Data input details
As any other simulation program (commercial or free), the PTL requires of information data as input, it is needed for the analysis process.The information can be given as a data file that contains basic information to generate the base study case: the base power of the system, nodal information (number of nodes, voltages, load powers and generation power), machine limits, system branches, SVC information (if applies).The simulation program PTL can handle two file extensions: cdf (standardized IEEE format) and ptl (proposed PTL format).

Simulator description
The input information for carrying out the search of the solution process, as the related data to the problem results of the load flow problem (for generating the base case) are stored in defined data structures (e.g.Dat_Vn, Dat_Gen, Dat_Lin, Dat_Xtr).These structures allow easily the data extraction, by naming each field in such a way that the programming becomes intuitive and each variable can be easily identified with the corresponding physical variable of the problem.An example with two structures used in the PTL program is shown in Table 1.The MATLAB structures are composed of non-primitive variable types that allow storing different data types in a hierarchical way with the same entity (García et al., 2005).They are formed by data containers called fields, which can be declared by defining the structure name and the desired field considering its value, e.g.Dat_Vn.Amp=1.02.

Output information
The PTL program displays the results obtained from the load flow execution in a boxlist (uicontrol MATLAB) with a defined format: nodal information, power flows in branches and generators.It gives the option of printing a report in Word format with the same information.In addition, it has the option of generating a file with the extension f2s which is oriented for stability studies and it contains all necessary information for the initialization of the state machine variables by using the results of the current power flow solution.
The definition of the conceptual simulation program PTL is presented as a recommendation by taking into account the three basic points that a simulation program must include: to be completely functional, to be general for any study case and to facilitate its maintenance; in other words it must allow the incorporation of new functions for the solution of new studies, such that it allows its free modification as easy as possible.

Load flows
In a practical problem, the knowledge of the operating conditions of an electric power system is always needed; that is, the knowledge of the nodal voltage levels in steady-state under loaded and generating conditions and the availability of its transmission elements are required to evaluate the system reliability.Many studies focused to the electric power systems start from the load flow solution which is known as "base case", and in some cases, these studies are used to initialize the state variables of dynamic elements of a network (generators, motors, SVC, etc) to carry out dynamic and transient stability studies.Another study of interest, that it also requires starting from a base case, is the analysis of the power system security that will be discussed in next sections.
The mathematical equations used to solve this problem are known as power flow equations, or network equations.In its more basic form, these equations are derived considering the transmission network with lumped parameters under lineal and balanced conditions, similarly as the known operating conditions in all nodes of the system (Arrillaga, 2001).

Power flow equations
An electric power system is formed with elements that can be represented for its equivalent circuit RLC, and with components as load and generating units which cannot be represented as basic elements of an electrical network, they are represented as nonlinear elements.However, the analysis of an electrical power system starts with the formulation of a referenced nodal system and it describes the relationship between the electrical variables (voltages and currents) as it is stated by the second Kirchhoff´s law or nodal law.
where I BUS is a n×1 vector whose components are the electrical net current injections in the n network nodes, V BUS is a n×1 vector with the nodal voltages measured with respect to the referenced node and Y BUS is the n× n nodal admittance matrix of the electrical network; it has the properties of being symmetric and squared, and it describes the network topology.
In a real power system, the injected currents to the network nodes are unknown; what it is commonly known is the net injected power S k .Conceptually, S k is the net complex power injected to the k-th node of the electrical network, and it is determined by the product of voltage (V k ) and the current conjugate (I k * ), where V k and I k are the voltages and nodal currents at the node k, that is, the k-th elements of vectors BUS V y BUS I in (1).Once the I k is calculated using (1), the net complex power S k can be expressed as: where , Y km is the element (k,m) of BUS Y matrix in (1).S k can also be represented for its real and imaginary components such as it is shown in the following expression: where P k and Q k are the net active and reactive power injected at node k of the system, respectively.They are defined as: where the variables  ( ) The equations ( 9) and ( 10) are commonly known as Power flow equations and they are needed for solving the load flow problem (Arrillaga, 2001).By analyzing these equations it can be clearly seen that each system node k is characterized for four variables: active power, reactive power, voltage magnitude and angle.Hence it is necessary to specify two of them and consider the remaining two as state variables to find with the solution of both equations.

Bus types in load flow studies
In an electrical power network, by considering its load flow equations, four variables are defined at each node, the active and reactive powers injected at node k P and k Q , and the magnitude and phase voltage at the node k V y k θ .The latter two variables determine the total electrical state of the network, then, the objective of the load flow problem consists in determining these variables at each node.The variables can be classified in controlled variables, that is, its values can be specified and state variables to be calculated with the solution of the load flow problem.The controlled or specified variables are determined by taking into account the node nature, i.e. in a generator node, the active power can be controlled by the turbine speed governor, and the voltage magnitude of the generator node can be controlled by the automatic voltage regulator (AVR).In a load node, the active and reactive power can be specified because its values can be obtained from load demand studies.Therefore, the system nodes can be classified as follows: • Generator node PV: It is any node where a generator is connected; the magnitude voltage and generated active power can be controlled or specified, while the voltage phase angle and the reactive power are the unknown state variables (Arrillaga, 2001).

•
Load node PQ: It is any node where a system load is connected; the active and reactive consumed powers are known or specified, while the voltage magnitude and its phase angle are the unknown state variables to be calculated (Arrillaga, 2001).

•
Slack node (Compensator): In a power system at least one of the nodes has to be selected and labeled with this node type.It is a generating node where it cannot be specified the generated active power as in the PV node, because the transmission losses are not known beforehand and thus it cannot be established the balance of active power of the loads and generators.Therefore this node compensates the unbalance between the active power between loads and generating units as specified in the PQ and PV nodes (Arrillaga, 2001).

Solution of the nonlinear equations by the Newton-Rhapson method
The nature of the load flow problem formulation requires the simultaneous solution of a set of nonlinear equations; therefore it is necessary to apply a numerical method that guarantees a unique solution.There are methods as the Gauss-Seidel and Newton-Rhapson, in the work presented here, the load flow problem is solved with the Newton-Raphson (NR) (Arrillga, 2001).
The NR method is robust and has a fast convergence to the solution.The method has been applied to the solution of nonlinear equations that can be defined as, where () fx is a n×1 vector that contains the n equations to be solved () () () , () , , () and x is a n×1 vector that contains the state variables, The numerical methods used to solve ( 11) are focused to determine a recursive formula The solution algorithm is based in the iterative application of the above formula starting from an initial estimate 0 x , until a convergence criterion is achieved ( ) , where ε is a small numerical value, and the vector k x is an approximation to the solution * x of ( 11).The Newton-Raphson method can be easily explained when it is applied to an equation of a single state variable (Arrillaga, 2001).A geometrical illustration of this problem is shown in figure 2. From figure 2, it can be established that where, As it can be seen in figure 2, the successive application of this correction 12 ,,. . .x as nearer as desired.To derive the recursive formula to be employed in the solution of the set of equations, and by expressing the equation ( 15) in matrix notation, it is obtained, where and the recursive formula is given by,

Application of the NR to the power flow problem
The equations to be solved in the load flow problem, as it was explained in section 3.1, are here rewritten, ( ) The 2n equations to be solved are represented by ( 4) and (5).However, for all generator nodes, equation ( 5) can be omitted and for the slack node the equations ( 4) and ( 5) (Arrillaga, 2001).The resulting set of equations is consistent because for the neglected equations, its corresponding state variables esp P , esp Q are also omitted from them.All this operation gives as a result a set of equations to solve, where its state variables x only contain magnitudes of nodal voltages and its corresponding phase angles which are denoted by V y θ, which simplifies considerably a guaranteed convergence of the numerical method.Therefore, the set of equation can be represented in vector form as, () where () ΔPx represents the equation ( 4) for PV and PQ nodes, () ΔQx represents the equation ( 5) for PQ nodes and x denotes the state variables V and θ, which are represented in vector notation as: where: Considering equation ( 16), its matrix equation is obtained and it defines the solution of the load flow problem: By applying the recursive equation ( 18), the state variables (V y θ) are updated every iteration until the convergence criterion is achieved max( ( ) ) ε Therefore, equation ( 21) can be expressed as, where the Jacobian elements are calculated using equations ( 9) and ( 10), provided equation ( 21) is normalized.
To simplify the calculation of the Jacobian elements of ( 24), it can be reformulated as: where the quotient of each element with its corresponding element V allows that some elements of the Jacobian matrix can be expressed similarly (Arrillaga, 2001).The Jacobian can be formed by defining four submatrixes denoted as H, N, J and L which are defined as, If k≠m (off diagonal elements): ( ) ( ) ( ) If k=m (main diagonal elements): It is important to consider that subscripts k and m are different, km km LH = and km km JN =− , on the contrary, to the principal diagonal elements.Nevertheless, there is a relationship when the following equation is solved: 11 00 00 00 00 00 00 By simple inspection of (34), we can see that the operation for km N − and km J results to be the real part of (34), while km H and km L are the imaginary component of it.The variant consists that all elements of the principal diagonal are calculated after (34) has been solved by subtracting or adding the corresponding k P or k Q as it can be seen in ( 30)-( 33).The Jacobian matrix has the characteristic of being sparse and squared with an order of the length of vector X, where its sub-matrixes have the following dimensions: H: n-1 x n-1 matrix.N: n-1 x nc matrix.J: nc x n-1 matrix.L: nc x nc matrix.

Voltage stability and collapse
The terms of voltage stability and collapse are closely related to each other in topics of operation and control of power systems.

Voltage stability
According to the IEEE, the voltage stability is defined as the capacity of a power system to maintain in all nodes acceptable voltage levels under normal conditions after a system disturbance for a given initial condition (Kundur, 1994).This definition gives us an idea of the robustness of a power system which is measured by its capability of keep the equilibrium between the demanded load and the generated power.The system can be in an unstable condition under a disturbance, increase of demanded load and changes in the topology of the network, causing an incontrollable voltage decrement (Kundur, 1994).The unstable condition can be originated for the operating limits of the power system components (Venkataramana, 2007), such as: Generators: They represent the supply of reactive power enough to keep the power system in stable conditions by keeping the standardized voltage levels of normal operation.However, the generation of machines is limited by its capability curve that gives the constraints of the reactive power output due to the field winding current limitation.Transmission lines: They are another important constraint to the voltage stability, and they also limit the maximum power that be transported and it is defined the thermal limits.Loads: They represent the third elements that have influence on the stability voltage; they are classified in two categories: static and dynamic loads and they have an effect on the voltage profiles under excessive reactive power generation.

Voltage collapse
The voltage collapse is a phenomenon that might be present in a highly loaded electric power system.This can be present in the form of event sequence together with the voltage instability that may lead to a blackout or to voltage levels below the operating limits for a significant part of the power system (Kundur, 1994).Due to the nonlinear nature of the electrical network, as the phenomenon related to the power system, it is necessary to employ nonlinear techniques for the analysis of the voltage collapse (Venkataramana, 2007) and find out a solution to avoid it.There many disturbances which contribute to the voltage collapse: To reach the reactive power limits in generators, synchronous condensers or SVC.

•
The operation of TAP changers in transformers.

•
The tripping of transmission lines, transformers and generators.Most of these changes have a significant effect in the production, consumption and transmission of reactive power.Because of this, it is suggested control actions by using compensator elements as capacitor banks, blocking of tap changers, new generation dispatch, secondary voltage regulation and load sectioning (Kundur, 1994).

Analysis methods for the voltage stability
Some of the tools used for the analysis of stability voltage are the methods based on dynamic analysis and those based in static analysis.Dynamic Analysis.They consist in the numerical solution (simulation) of the set of differential and algebraic equations that model the power system (Kundur, 1994), this is similar as transients; however, this kind of simulations need considerable amount of computing resources and hence the solution time is large and they do not give information about the sensibility and stability degree.Static analysis.They consist in the solution of the set of algebraic equations that represent the system in steady state (Kundur, 1994), with the aim of evaluating the feasibility of the equilibrium point represented by the operating conditions of the system and to find the critical voltage value.The advantage with respect to the dynamic analysis techniques is that it gives valuable information about the nature of the problem and helps to identify the key factors for the instability problem.The plotting of the PV curve helps to the analysis of the voltage stability limits of a power system under a scenario with load increments and with the presence of a disturbance such as the loss of generation or the loss of a transmission line.

PV curves
The PV curves represent the voltage variation with respect to the variation of load reactive power.This curve is produced by a series of load flow solutions for different load levels uniformly distributed, by keeping constant the power factor.The generated active power is proportionally incremented to the generator rating or to the participating factors which are defined by the user.The P and Q components of each load can or cannot be dependant of the bus voltage accordingly to the load model selected.The determination of the critical point for a given load increment is very important because it can lead to the voltage collapse of the system.These characteristics are illustrated in figure 3. Some authors (Yamura et al, 1998& Ogrady et al, 1999) have proposed voltage stability indexes which are based in some kind of analysis of load flows with the aim to evaluate the stability voltage limits.However, the Jacobian used in the load flows, when the Newton-Raphson is employed, becomes singular at the critical point, besides the load flow solutions at the points near to the critical region tend to diverge (Kundur, 1994).These disadvantages are avoided by using the method of continuation load flows (Venkataramana et al, 1992).

Application of the continuation method to the power flow problem
The continuous load flow procedure is based in a reformulation of the equations of the load flow problem and the application of the continuation technique with a local parameterization which has shown to be efficient in the trajectory plotting of PV curves.The purpose of continuous load flows is to find a set of load flow solutions in a scenario where the load is continuously changing, starting from a base case until the critical point.Thereafter, the continuous load flows had been applied to understand and evaluate the problem of voltage stability and those areas that are likely to the voltage collapse.Besides, they have also been applied in other related problems like the evaluation of power transfer limits between regions.The general principle of continuous load flows employs a predictor-corrector scheme to find a trajectory of solutions for the set of load flow equations ( 4) and ( 5) which are reformulated to include the load parameter λ.

(
) The process is started from a known solution and a predictor vector which is tangent to the corrected solutions is used to estimate the future solutions with different values of the load parameter.The estimation is corrected using the same technique of the Newton-Rhapson employed in the conventional load flow with a new added parameter: ( ) The parameterization plays an important role in the elimination of the Jacobian nonsingularity.

Prediction of the new solution
Once the base solution has been found for λ=0, it is required to predict the next solution taking into consideration the appropriate step size and the direction of the tangent to the trajectory solution.The first task in this process consists in calculating the tangent vector, which is determined taking the first derivative of the reformulated flow equations ( 38).
(, , ) 0 where F is the vector [ΔP, ΔQ, 0] that is augmented in one row; in a factorized form, the equation is expressed as, [ ] The left hand side of the equation is a matrix of partial derivatives that multiplies the tangent vector form with differential elements.The matrix of partial derivatives is known as the Jacobian of the conventional load flow problem that is augmented by the column F λ , which can be obtained by taking the partial derivatives with respect to λ (35) and ( 36), which gives: Due to the nature of (41) which is a set of nc+n-1 equations with nc+n unknowns and by adding λ to the load flow equations, it is not possible to find a unique and nontrivial solution of the tangent vector; consequently an additional equation is needed.This problem is solved by selecting a magnitude different from zero for one of the components of the tangent vector.In other words, if the tangent vector is denoted by: where e k is a vector of dimension m+1 with all elements equal to zero but the k-th one, which is equal to 1.If the index k is correctly selected, t k = ±1 impose a none zero norm to the tangent vector and it guarantees that the augmented Jacobian will be nonsingular at the critical point (Eheinboldt et al, 1986).The usage of +1 or -1 depends on how the k-th state variable is changing during the trajectory solution which is being plotted.In a next section of this chapter, a method to select k will be presented.Once the tangent vector has found the solution of ( 43), the prediction is carried out as follows: where "*" denotes the prediction for a future value of λ (the load parameter) and σ is a scalar that defines the step size.One inconvenient in the control of the step size consists in its dependence of the normalized tangent vector (Alves et al, 2000), where |t k | is the Euclidian norm of the tangent vector and σ 0 is a predefined scalar.The process efficiency depends on making a good selection of σ 0 , which its value is system dependant.

Parameterization and corrector
After the prediction has been made, it is necessary to correct the approximate solution.
Every continuation technique has a particular parameterization that gives a way to identify the solution along the plotting trajectory.The scheme here presented is referred as a local parameterization.In this scheme, the original set of equations is augmented with one extra equation, which has a meaning of specifying the value of a single state variable.In the case of the reformulated equations, this has a meaning of giving a unity magnitude to each nodal voltage, the phase angle of the nodal voltage or the load parameter λ.The new set of equations involves the new definition of state variables as: where, where k x ε x and η represent the appropriate k-th element of x.Therefore, the new set of equations that substitute (38) is given by, () 0 After an appropriate index k has been selected and its corresponding value of η , the load flow is solved with the slightly modified Newton Raphson method (48).In other words, the k index used in the corrector is the same as that used in the predictor and η is equal to the obtained x k from the corrector (44), thus the variable x k is the continuation parameter.
The application of the Newton-Raphson to (38) results in, where e k is the same vector used in (43), and the elements ΔP and ΔQ are calculated from ( 35) and (36).Once the x k is specified in ( 48), the values of the other variables are dependant on it and they are solved by the iterative application of (49).

Selection of the continuation parameter
The most appropriate selection corresponds to that state variable with the component of the tangent vector with the largest rate of change relative to the given solution.Typically, the load parameter is the best starting selection, i.e. λ=1.This is true if the starting base case is characterized for a light or normal load; in such conditions, the magnitudes of the nodal voltages and angles keep almost constants with load changes.On the other hand, when the load parameter has been increased for a given number of continuation steps, the solution trajectory is approximated to the critical point, and the voltage magnitudes and angles probably will have more significant changes.At this point, λ has had a poor selection as compared with other state variables.Then, once the first step selection has been made, the following verification must be made: where m is equal to the state variables; including the load parameter and k corresponding to the maximum t/x component.When the continuation parameter is selected, the sign of the corresponding component of the tangent vector must be taken into account to assign +1 or -1 to t k in (42) for the subsequent calculation of the tangent vector.

MATLAB resources for electrical networks
The reason why MATLAB is frequently chosen for the development of academic or research tools is because its huge amount of mathematical operations as those related to vectors and matrixes.In addition, it also has several specialized libraries (toolboxes) for more specific areas as control, optimization, symbolic mathematics, etc.In the area of power systems, it is possible to point it out several advantages considered as key points for the development of a script program, which are discussed below.

Sparse matrix manipulation
The electrical networks are studied using nodal analysis, as it was presented in section 3.1.
In this frame of reference, the network matrixes as Ybus or Jacobian (1) and ( 21) have a sparse structure, considering that a node in an electrical network is connected to about 2.4 nodes in average.To illustrate this issue, if we consider the creation of a squared matrix with order 20, i.e. the matrix has 200 elements which 256 are zeros and the remaining are different from zero that are denoted as nz, as it is shown in figure 4. The sparsity pattern is displayed by using the MATLAB function spy.The operation on this type of matrices with conventional computational methods leads to obtain prohibited calculation times (Gracia et al, 2005).Due to this reason, there have been adopted special techniques for deal with this type of matrixes to avoid the unnecessary usage of memory and to execute the calculation processes on the nonzero elements.It is important to point it out that this kind of matrixes is more related to a computation technique than to a mathematic concept.MATLAB offers mechanisms and methods to create, manipulate and operate on this kind of matrixes.The sparse matrixes are created with the sparse function, which requires the specification of 5 arguments in the following order: 3 arrays to specify the position i, the position j -row and column-and the element values x that correspond to each position (i,j), and the two integer variables to determine the dimensions mxn of the matrix, for example: >> A = sparse(i, j, x, m, n) In MATLAB language, the indexing of full matrixes is equal to the sparse matrixes.This mechanism consists in pointing to a set of matrix elements through the use of two arrays that makes reference to each row and column of the matrix, e.g.B = A[I, J] or by simply pointing to each element o elements of the matrix to be modified, e.g.A[I, J] = X.
where A is a mxn spare matrix, I and J are the arrays that point to the rows and columns and X is the array that contains each element corresponding to each ordered pair (I k , J k ).All mentioned arrays are composed for a number of k elements, such as k < n and k < m.
MATLAB has a function called spdiags that is intended for the direct operation on any diagonal of the matrix, e.g. to make uniform changes on the elements of the diagonal "d" of matrix A: >> A = spdiags(B, d, A) A mechanism frequently used to form matrixes from other defined ones is called the concatenation.Even more, to form sparse matrixes by using other arrays of the same type but with specified dimensions in such a way that there is a consistency to gather into a single one, i.e. the concatenating is horizontal, the matrix rows must be equal to those of B matrix, e.g.C = [A B].On the other hand, if the concatenating is vertical, then the columns of both matrices must be equal, e.g.C = [A; B].

The LU factorization for solving a set of equations
In electrical network applications, to find a solution of the algebraic system = Ax b in an efficient way, one option is to use the triangular decomposition as the LU factorization technique.The triangular factorization LU consists in decomposing a matrix (A) such that it can be represented as the product of two matrixes, one of them is a lower triangular (L) while the other is an upper triangular (U), =⋅ AL U .This representation is commonly named explicit factorization LU; even though it is very related to the Gaussian elimination, the elements of L and U are directly calculated from the A elements, the principal advantage with respect to the Gaussian elimination (Gill et al., 1991) consists in obtaining the solution of an algebraic system for any b vector, if and only if the "A" matrix is not modified.

Sparse matrix ordering using AMD
The Approximate minimum degree permutation (AMD) is a set of subroutines for row and column permutation of a sparse matrix before executing the Cholesky factorization or for the LU factorization with a diagonal pivoting (Tim, 2004).The employment of this subroutine is made by using the function "amd" to the matrix to be permuted, e.g.P = amd(A).

Sparse matrix operations
Care must be taken with a several rules in MATLAB when operations are carried out that include full and sparse matrixes, for example, the operation eye(22) + speye(22) gives a full matrix.

Vector operations
An approach to develop the script programs in MATLAB to be executed faster consists in coding the algorithms with the use of vectors within the programs and avoiding the use of loops such as for, while and do-while.The vector operations are made by writing the symbol "." before the operation to be made, e.g..+,.-,.*,./.This discussion is illustrated by comparing two MATLAB script programs to solve the following operation: power transfer, etc. Therefore the input data related to these controllers is not used in the program.

Graphical interface for editing the display of PV curves
In the last part of the program, a graphical user interface (GUI) is generated for the plotting of PV curves from the points previously calculated in the continuous load flows program.
Even though the main objective of this interface is oriented to the plotting of PV curves, its design allows the illustration continuation method process, i.e. it has the option of displaying the calculated points with the predictor, corrector or both of them, for example, it can be seen in figure 5 the set of points obtained with the predictor-corrector for the nodes 10,12,13, and 14 of the 14-node IEEE test power system.This is possible by integrating the MATLAB controls called uicontrol´s into the GUI design, the resultant GUI makes more intuitive and flexible its use.
The numbers 1-12, that are indicated in red color in figure 5, make reference to each uicontrol or graphical object that integrates the GUI.The graphical objects with its uicontrol number, its handle and its description are given in table 1.Since the function PV_PRINT is external to the function main, the created variables within main cannot be read directly for the function PV_PRINT, hence it is required to declare arguments as inputs for making reference to any variable in main.The main advantage of using external functions consists in the ability to call them in any application just by giving the required arguments for its correct performance.In order to reproduce this application, it is necessary to copy the code in an m-file which is different from the main.

Conclusion
An electrical engineering problem involves the solution of a series of formulations and mathematical algorithm definitions that describe the problem physics.The problems related to the control, operation and diagnostic of power systems as the steady-state security evaluation for the example the PV curves, are formulated in matrix form, which involves manipulation techniques and matrix operations; however the necessity of operating on matrixes with large dimensions takes us to look for computational tools for handling efficiently these large matrixes.The use of script programming which is oriented to scientific computing is currently widely used in the academic and research areas.By taking advantage of its mathematical features which are normally found in many science or engineering problems allows us solving any numerical problem.It can be adapted for the development of simulation programs and for illustrating the whole process in finding a solution to a defined problem, and thus makes easier to grasp the solution method, such as the conventional load flow problem solved with the Newton-Raphson method.MATLAB has demonstrated to be a good tool for the numerical experimentation and for the study of engineering problems; it provides a set of functions that make simple and straightforward the programming.It also offers mechanisms that allow dealing with mathematic abstractions such as matrixes in such a way that is possible to develop prototype programs which are oriented to the solution methods by matrix computations.The development of scripts or tools can be considered to be a priority in the academic area such that they allow achieving a valid solution.It is also advisable to take advantage of the MATLAB resources such as: vector operations, functions and mechanisms for operating on each matrix element without using any flow control for the program, i.e. for loop; this offers the advantage of decreasing the number of code lines of the script program.These recommendations reduce the computation time and allow its easy usage and modification by any user.Finally, MATLAB offers powerful graphical tools which are extremely useful for displaying the output information and to aid interpreting the simulation results.In this chapter, the plotting of the corrector points (PV curve) has been presented for a given load

Fig. 5 .
Fig. 5. GUI for plotting the PV curves for any node.

Table 1 .
Example of the information handled in the PTL.