TESTING A DIFFERENTIAL-ALGEBRAIC EQUATION SOLVER IN LONG-TERM VOLTAGE STABILITY SIMULATION

This work evaluates the performance of a particular differential-algebraic equation solver, referred to as DASSL, in power system voltage stability computer applications. The solver is tested for a time domain long-term voltage stability scenario, including transient disturbances, using a real power system model. Important insights into the mechanisms of the DASSL solver are obtained through the use of this real model, including control devices relevant to the simulated phenomena. The results indicate that if properly used, the solver can be a powerful numerical tool in time domain assessment of long-term power system stability since it comprises, among several important features, suitable and very efficient variable order and variable step-size numerical techniques. These characteristics are very important when CPU time is a great concern, which is the case when the power system operator needs reliable results in a short period of time. Prior to the present work, this solver has never been applied in power system stability computer analysis in time domain considering slow and fast phenomena.


Introduction
Time domain techniques for power system stability studies have been kept under investigation due to the difficulties that arise when one links complex mathematical models and efficient numerical methods.Power system time domain computer simulation requires, in general, the solution of large-scale differential-algebraic equation systems, DAEs, which, besides the complexity of solving these DAEs systems, may be computationally time consuming in terms of CPU.This computational effort depends mainly on the available numerical algorithm, complexity of the mathematical models, power system dimension, dynamic-devices time-constant range, speed of the phenomenon (fast or slow), computer capability, and total time simulated which may range from few seconds 2 Long-term voltage stability simulation-DAEs solution up to tens of thousands of seconds.The present work deals with two specific forms of power system stability, named as synchronous machine angular transient stability and long-term voltage stability.
Synchronous machine angular transient stability phenomenon is, in general, numerically solved by implicit fixed-step-size multistep methods, such as trapezoidal integration with Newton iteration.Depending on the size of the system being solved, these methods can be time consuming because they do not lead to accurate approximations for large integration step sizes.The step size must remain within the region of absolute stability throughout the interval of t values for the approximation decay to zero and the growth error to be under control [4].
On the other hand, time domain simulation of long-term phenomena using fixed and small step sizes is computationally prohibitive.Slow acting control devices, such as overexcitation limiters and/or under load tape changers, can take tens/hundreds of seconds to act.Variable step size methods are especially suitable in these cases, since an examination of the local truncation error might indicate that the step size could increase.Despite of the actual and efficient variable step size algorithms [5], time domain power system stability analysis still represents a challenging task in computer modeling and simulation for real time operation.
In time domain power system studies, the use of simplified models may reduce computer time processing, however, depending on the scenario being simulated, the accuracy may be lost since relevant effects are not captured.The analysis can still become more complex if different forms of stability come about simultaneously.On the other hand, very complex models demand a considerable effort in the solution of the DAEs and this complexity may be unnecessary for some cases.Therefore, the use of suitable DAEs techniques in long-term power system voltage stability studies is a very important issue deserving special attention.Some DAEs techniques have been proposed and tested in different scientific fields and some of them are implemented as computer codes, known as solvers.The basic differences among them are the numerical tools offered by each solver.In function of that, some solvers present better performance for specific applications, and for time domain power system stability studies, not much has been found in the literature [1,6,12].
This work deals with a differential-algebraic system solver, known as DASSL [7], which is based on BDF (backward differentiation formulae) methods, variable step size, and variable order.The solver has been designed to be used for the solution of DAEs of index zero and one and uses the fixed-leading-coefficient form of the BDF formulae to solve general index-1 DAEs.The DASSL Fortran version used in this work is a public domain software [2,3] (http://www.engineering.ucsb.edu/%7Ecse/).

Origin of the problem
The basic structure representing the power system for voltage stability analysis in time domain is similar to that used in conventional analysis of angular transient stability, where mathematical models are represented by differential-algebraic equations.Power system voltage stability time domain analysis is useful in determining the time coordination of the equipments; clarifying the phenomenon and preventing overdesign, improving simulation fidelity, time domain simulation forces a more careful analysis and the use (development) of more perfect models; simulating fast dynamics associated with the final phases of the collapse; and obtaining the dynamic acting of the system through graphs where the evolution is visualized in the time of the phenomenon of the angular/voltage stability [13].
Power system angular transient stability in time domain is often analyzed based on the behavior of the synchronous machine rotor's angular position (δ) related to a reference axis running in synchronous speed, between two or more synchronous machines.Although the angular position corresponds to the variable of interest for transient angular stability evaluation in time domain, important information related to transient voltage stability (fast voltage phenomena) can be obtained through the mathematical models that describe the dynamic behavior of synchronous machines.Conventional transient stability (angular) programs can be used for transient voltage stability analysis, certifying that models relevant to the phenomenon are available, such as induction motors and directcurrent converters.Voltage collapse can occur as a consequence of a transient voltage instability condition due to unfavorable response of these devices to a system disturbance [13].
However, care must be taken when traditional methods for transient stability analysis are used in long-term voltage stability studies in time domain.The generator's overexcitation limiter (OXL) is an important control device for long-term voltage stability studies that is not considered in transient stability studies.The effect of the OXL over the generator's terminal voltage is important and this device acts in a period of time higher than the transient period (t ≥ 10 s).Therefore, this device is not considered or other devices of slower dynamics, in conventional studies of angular stability in time domain.
The period of time concerned in the study can be extended and the efficiency and the limits of this extension depend mainly on the numerical method used in the solution of the DAEs and on the models representing control devices, loads, and other equipment.If the study is extended to the minutes range, the slower dynamic devices, if available, can detect unfavorable conditions.This is a very important issue in long-term voltage stability studies.
Voltage stability problems can lead the power system to a partial or total voltage collapse.This scenario can last tens of minutes (long-term) but fast dynamics should not be neglected.Transients phenomena can occur as a result of severe disturbances, such as losses of generation units or openings of important transmission lines.Regarding this statement, it is important to note the following.
(i) Fast dynamics are not important only at the beginning (initial disturbance) or at the end (voltage collapse final form) of the simulation.(ii) The number of transients occurrences is related to the number of severe disturbances, probably occurring several times in the same simulation (sequence of events).
The above statements clarify the necessity of a reliable and efficient variable-step-size algorithm if slow and fast phenomena are to be simulated, since they range from few seconds up to tens of thousands of seconds.

The DAE solver
There are several DAEs solvers and some of them are freely available computer codes (test system for IVP ODEs, http://www.unige.ch/math/folks/hairer/).The great advantage of these solvers is that, in general, they are ready to use and comprise efficient numerical techniques.Besides, if the computational code source is made available, the user can introduce suitable changes adjusting the solver according to his needs.In this work, the DASSL code was chosen because, at first view, it presented important and suitable numerical features that seemed very adequate for long-term voltage stability applications.Some of these features are: (1) variable-order and variable-step-size algorithm; (2) the method implemented in (3.1) belongs to the BDF class, having great ability in solving stiff differential-algebraic equations [2].The DAEs system representing the power system dynamic devices including fast and slow response control devices is stiff; (3) availability of a good literature and computational code source; (4) originality for long-term voltage stability applications.Prior to the present work, this solver has never been applied in power system long-term voltage stability time domain analysis.Regarding the first issue above, varying the step size is essential for the effective performance of a discretization method.The problem is how the DAE code decides to accept or reject a step.In general, the errors made at each step are much easier to estimate than the global error.Even though the global error is more significant, the local truncation error is the one that general-purpose multistep codes usually estimate in order to control the step size and to decide on the order of the method to be used.DASSL deals with two sources of errors.The first source of error treated by DASSL is the local truncation error of the method, estimated directly by using divided differences.The second source is the error in interpolating to find the solution between mesh points for output purposes [3].
DASSL applies a variable-step-size variable-order fixed-leading-coefficient implementation of BDF formulas to advance the solution from one time step to the next.The fixed leading coefficient is an interesting strategy of extending the fixed-step-size BDF methods to variable step sizes.The fixed leading coefficient includes the best features of fixedcoefficient strategy and variable-coefficient strategy.The fixed-coefficient strategy has the advantage of simplicity, but they are less stable than variable-coefficient formulae due to an error inherent in the interpolation.Stability of the variable step size formulae is an essential concern for problems where the step size must be changed frequently or drastically [3].
The variable-coefficient methods are the most stable implementation, but they tend to require more evaluations of the Jacobian matrix in intervals when the step size changes, and therefore are usually considered to be less efficient than a fixed coefficient implementation for most problems.This method has the advantage for problems which call for repeated or radical changes of the step size.This is an appropriate characteristic for time domain analysis including in the same simulation fast and slow phenomena.
Another important attribute is that DASSL possesses an interpolant to compute the solution at points different from time steps chosen by the algorithm.This is an important Figure 3.1.Brazilian south-southeast power system equivalent model.aspect for power system stability simulation since the time sequence concerning the events are reproduced with a good fidelity.The loss of accuracy may result, for instance, in incorrect information for relays settings and other control devices adjustments.This may lead one or more protective devices to an inappropriate operation, triggering a detrimental sequence of events, whose final result may be catastrophic (blackouts).

Testing the solver.
The solver performance is tested for transient and long-term angular/voltage stability phenomena using the power system model illustrated in Figure 3.1.This power system represents the Brazilian south-southeast equivalent system (private communication with Brazilian Electrical Energy Research Center, CEPEL) and comprises forty-two (42) buses.Reactive power shunt compensation devices are not represented.The main reason in testing the code for this kind of problem is based on the fact that simulating short-and long-term power system dynamics is not a straightforward computational task.In this particular case, the resulting DAEs system representing the network dynamic models is stiff.Besides, the numerical algorithm must recognize the right moment for reducing the integration step, as well as the moment for increasing it.In response to a series of power system internal contingencies, the dynamic trajectory of the operating point may present a large number of discontinuities.In real life, depending on specific characteristics, such as local and amplitude of the disturbance, fast and/or slow control devices responses, and operator actions, the power system may take few seconds, or even hours to reach, or not, a stable operating condition.The power system may return 6 Long-term voltage stability simulation-DAEs solution to the initial operating condition, or reach a different one.Therefore, a variable step-size and robust algorithm for solving stiff DAEs systems is the right choice and the DASSL code seems much appropriated in this case.

Mathematical description of the problem.
The problem is a stiff DAE of index 1, consisting of 702 equations having the following form: with y ∈ 702 , 0 ≤ t ≤ 2783 (t in seconds).
For this specific problem, there are 548 differential equations and 154 algebraic equations, and M is the zero matrix except for M 1•••548;1•••548 , which is given by To illustrate the differential and algebraic equations placed into the solver, we present bellow the group of equations representing the synchronous machine dynamic model and respective control devices, and the algebraic equations representing the network.The notation is standard for a synchronous machine model including transient and subtransient effects, and saturation as well.It is assumed that the stator/network electromagnetic transients have been eliminated, leading to algebraic equations that accompany the multimachine dynamic model.Besides, the stator/network electromagnetic transients are very fast in comparison with the angular and voltage phenomena under investigation.They can be eliminated using singular perturbation and the concept of integral manifolds [8][9][10].
(i) mth synchronous machine internal equations: (ii) mth synchronous machine stator algebraic equation: (iii) mth synchronous machine inductances: (iv) mth synchronous machine automatic voltage regulator output (AVR-Efd): (v) mth synchronous machine power system stabilizer output (PSS-V out ): 8 Long-term voltage stability simulation-DAEs solution (vi) mth synchronous machine overexcitation limiter output (OXL-V OXL ): (vii) Network algebraic equations (including generation buses m and load buses n): Taking into account the above differential and algebraic equations, the corresponding DAE system representing the electrical power system and its components is formed by [11] (a) five differential equations for each machine: (b) two stator algebraic equations of the form (3.10) in the polar form:  The simulation consists of a line opening between buses 28 and 29 at t = 5 seconds, followed by a linear reactive power (MVAr) load pick up at a rate of 0.035 pu s −1 in bus 29.Overexcitation limiters of all generators have been adjusted to operate for a field current of 10% above their respective nominal values.Remembering that all synchronous generators comprise AVR and PSS, which are fast response control devices, and together with generators' OXL and power system's ULTCs, which are slow response control devices, makes possible the simulation of transient and long-term phenomena.
As the Mvar load increases, the generators' field currents also increase, as shown in Figure 3.2.The G26 OXL reaches its set point first bringing the field current closer to its nominal value.The remaining generators assume this extra Mvar demand for the next seconds until the OXL of each synchronous generators group act.All generators' terminal voltages are no longer kept constant by the AVRs and the voltage level is further deteriorated, as can be seen in Figure 3.3.In addition to the OXL operation, the ULTCs act in order to keep constant the voltage level in buses 27 and 30, as shown in Figures 3.4 and 3.5.Not much can be done to control voltage levels when these devices reach their limits.The step size behavior is depicted in Figure 3.7.The discontinuities, characterized by the drastic step size reduction, are caused by the OXLs and ULTCs.This illustrates the solver capability of dealing with discontinuities caused by control devices.

Conclusions
This paper tested a free domain software for solving index zero and one systems of differential-algebraic equations, named as DASSL.The code encompasses an efficient variable step size and variable order based on BDF methods to solve a system of DAEs or ODEs.The computer code was applied in time domain power system stability studies.Using a real power system model including fast and slow control devices, it was possible to investigate the code capability of simulating different stability phenomena in the same run.The variablestep-size and variable-order algorithm implemented in DASSL is a very powerful tool for power system time domain computer simulation.The effects related to synchronous machines angular transient stability as well as to long-term voltage stability were simulated and captured in an efficient and suitable manner.Besides, the step sizes reach a range of tens of seconds, reducing the processing time.
There are advantages of applying DASSL in time domain studies, for instance, the code is ready to use, there is a public domain version available in the authors' code homepage (http://www.engineering.ucsb.edu/%7Ecse/), the solver background and techniques are very well documented in the literature, and it is applicable in a variety of problems (not only in power system time domain studies).

Table 3 .
1 presents the run characteristics, performed on a Pentium IV PC, 1.6 GHz, using Visual Fortran with optimization where, NSTEP denotes the number of time steps J. E. O. Pessanha and A. A. Paz 11 , NJAC the number of Jacobian evaluations, and CPU time the total CPU time taken to solve the problem, in seconds.From this table, it can be noted that the CPU time is very small, even taking into account the system's size.Several dynamic devices with different time frames are being considered, with time constant ranging from milliseconds up to 12 Long-term voltage stability simulation-DAEs solution * Figure 3.5.Controlled voltage at bus 30.used