Analysis of Drawbacks and Constraints of Classification Algorithms for Three-Phase Voltage Dips

Voltage events are one of the most common and harmful disturbances of power electric systems. Voltage dips, swells and interruptions are included under this heading. Given the economic cost that these disturbances represent for electrical power transmission and distribution companies and the industry, it becomes imperative to detect and classify them properly. Several classification criteria and algorithms have been proposed in the literature as analysis tools to differentiate voltage events by their characteristics and, if possible, to determine their causes and consequences. Even though some of these approaches make a correct classification of the voltage events, there are certain operation conditions that are common in real electrical grids, in which the classification criteria, and their corresponding algorithms, make a wrong classification. These particular conditions, together with the lack of a fair comparison in a common scenario, have not been addressed in the specific field literature. This work explores in detail all these aspects by evaluating the symmetrical components criterion and ABC classification criterion, and rigorously analyzes three specific algorithms: the Symmetrical Components Algorithm, the Six Phases Algorithm and the Space Vector Algorithm. Drawbacks arise from both classification criteria and algorithms. The causes of the classification errors are described and discussed in detail in order to better understand the problem, and evidence the constraints of these classification methods. Copyright © 2018 Praise Worthy Prize S.r.l. All rights reserved.

Nomenclature α = arg(Z f /Z s )
Phase angle between impedances β Argument of the V β phasor V a1 ,V b1 ,V c1 Phasors of the three-phase voltages φ +1 Phase angle of the positive sequence component of the fundamental voltage φ −1 Phase angle of the negative sequence component of the fundamental voltage φ inc Inclination angle of the ellipse formed by trajectory of the space vector a Unit phasor whose phase angle is 120° A, B, C, …,I * * Three-phase voltage dips of the ABC classification b Unit phasor whose phase angle is 60° D a , D b , D c Single-phase voltage dips of the CS classification C a , C b , C c Phasor whose argument β allows to estimates the T parameter Introduction In the Power Quality field, voltage disturbances are defined as deviations of the actual voltage signal from its ideal counterpart, where the ideal voltage is a sinusoidal waveform with constant amplitude and frequency parameters equal to the nominal ones [1]. One of the most harmful disturbances are events, which are characterized by a perfectly defined start and end times, such as dips, swells and interruptions. According to IEEE 1159-2009 Standard [2], voltage dips are defined as a reduction of the RMS (Root Mean Square) voltage value from 0.1 up to 0.9p.u. (per unit of the nominal or sliding reference voltage) in any of the three phases, with a duration from half-cycle up to one minute. This standard also defines a swell as the increment in the RMS value of the voltage above 1.1 p.u. in any of the three phases with a duration from half-cycle up to one minute; and a voltage interruption as a decrease of the RMS value of the voltage below 0.1 p.u. during less than a minute in any of the three-phases. Even though voltage dips are not harmless as interruptions or swells, they occur more often. Then, the damage produced by the voltage reduction, mainly given by short-circuits, motor starting and transformer energizing, can be higher [3]. As a consequence, their detection and classification have become a key issue for power systems operation, since knowledge about these events allows to determine the type of fault [4], and the source of disturbance and its location [5], among others, and to take corrective actions by clearing the fault, mitigating its effects [6] or avoiding repetition. The detection and classification of a voltage event can be automatically performed, by adopting a suitable classification criterion and implementing the corresponding algorithm in a digital signal processor.
The classification criteria allow determining the type of voltage dip, measuring its severity and providing information about the possible causes. Among the most popular classification criteria are the symmetrical components (SC) and the ABC criteria [7]. SC criterion is more general, since it classifies unbalanced dips as single, double or three-phase dips. Then, the dips produced in an actual network mainly fall within either group. On the other hand, the ABC criterion is more descriptive and includes a wide variety of events, based on more specific morphological features. For this reason, the ABC classification criterion is the most widely used, and is adopted in this manuscript. Nonetheless, both criteria are related, as indicated throughout this work.
Once the classification criterion is defined, the distinctive characteristics of the event that permit its identification should be extracted with an algorithm.
Some state-of-the-art works contribute with a theoretical approach about classification algorithms based on heuristic tools, such as fuzzy logics, neural networks, support vector machines, and genetic algorithms, to name a few [8]- [12]. Even though these methods have shown good performance in computer simulations, few experimental verification is provided in the cited works.
In addition, the implementation of these methods requires some particular skills like, for example, a good selection of the model and training sets, which could render them unattractive for industrial applications in which a known and simple method is most preferred. On the other hand, a widely adopted technique is the use of logical algorithms since they can be easily implemented with "nested if" structures.
This key feature justifies the analysis of the following event classification algorithms: the symmetrical components algorithm (SCA) [13], the six-phase algorithm (SPA) [13] and the space vector algorithm (SVA) [14]. Both SCA and SPA were designed to classify dips according to SC classification criterion, whose performance is described in [13].
Despite the fact that some suggestions for classification of voltage dips using ABC criterion from SC criterion and estimation of additional parameters are introduced in [15], the algorithms are not implemented and their performance and drawbacks are not analyzed.
On the other hand, the space vector algorithm is designed for ABC classification. In this case, the drawbacks of the algorithm and its performance under an actual event scenario is not available in the literature [14], [16]- [20]. As regards the algorithm operation, it is demonstrated in the literature that two of the three algorithms evaluated in this work (SCA and SPA) may fail under certain conditions [13], but no comments are made with respect to the third algorithm (SVA). Despite the fact that errors under particular conditions are reported in some works, they are not fully analyzed. In addition, while a topic such as power quality is widely discussed, the specific literature related to this kind of classification criteria and algorithms remains scarce.
In this work, the distinctive characteristics of the two criteria and the three classification algorithms are identified, the concepts related to classification criteria and algorithms are unified, and a fair comparison under the same operational conditions is conducted.
Furthermore, the characteristics and errors related to each criterion are thoroughly discussed, in order to evidence the real causes of classification errors. It is demonstrated that misclassifications result from not only the algorithms themselves, but also from the classification criterion on which they are based.
Moreover, the errors introduced when the voltage reference is different from the nominal voltage are analyzed. Finally, some aspects related to the possible improvement of the performance of these algorithms are discussed.
Section II introduces the main classification criteria. Then, the three algorithms under study are presented in Section III, while Section IV describes the limitations of both the algorithms and classification criteria, together with the causes of such limitations. Based on the concepts previously introduced, the algorithms are evaluated with simulations in Section V. Finally, the issues related to the causes leading to criteria and algorithm misclassification are dealt with in Section VI, after which conclusions are drawn in Section VII.

II. Classification Criteria of Voltage Dips
A common way of characterizing the severity of a three-phase voltage dip is to determine the lowest RMS value of the three phases during the event, and the highest time interval in which at least one phase is below the threshold of 0.9p.u. [21]. However, this methodology presents some limitations [22]. It classifies single, double or three-phase voltage dips with the same severity, which makes it difficult to draw conclusions regarding the type and location of the fault and does not provide clear information to analyze how a voltage dip propagates down to the power grid, for instance when passing through a transformer. To overcome these limitations, two classification criteria that provide a suitable characterization of the voltage dip are introduced in the work by Bollen and Zhang [7], known as Symmetrical Components (SC) and ABC.

II.1. Symmetrical Components Classification Criterion
The Symmetrical Components (SC) classification criterion, based on the analysis of voltage symmetrical components, distinguishes between voltage dips with the main voltage decrease in one phase (type D) and voltage dips with the main voltage decrease in two phases (type C), as shown in Figs. 1. Dips type D can be classified into Da, Db or Dc, where the subscript indicates the most affected phase. On the other hand, dips type C can be classified into Ca, Cb or Cc, where the subscript indicates the less affected phase. According to this classification, voltage dips can be represented as a phasor with two complex parameters, V and F. Factor V , known as characteristic voltage, depends on the impedances between the PCC (Point of Common Coupling), the fault location and the sources (see Eq. (1-A) in Appendix). Then, factor V should coincide with the lowest RMS voltage, which could be a line-to-line voltage affected by the factor √3 as shown in Fig. 1(a); or a single voltage, as shown in Fig. 1(b). The complex factor F, known as the PN (Positive Negative) factor, has a module that depends on the difference between the positive and negative sequence impedances located between the source and the PCC [7]. Each type of dip is defined by a set of equations expressed as a function of V and F. For instance, the definition of a type Ca is:

II.2. ABC Classification Criterion
The ABC classification criterion presented in [1] distinguishes seven types of dips (A, ..., G), and helps to link the event to the originating fault. The method was developed as part of a stochastic prediction of voltage dips. When a statistical analysis of faults is performed, it is possible to calculate how often events occur.
Moreover, the ABC classification allows to synthetize a wide variety of voltage dips; this is particularly useful for testing three-phase equipment that could operate under any of these events when connected to the electrical network. In [23], this criterion has been proposed as a complement to the tests recommended in the IEC 61000-4-11 Standard [24].
Originally, this classification considered a symmetrical voltage dip, A, three types of dips whose main voltage decrease was in one phase, B, D and F , and three kinds of dips whose main voltage reduction was in two phases, C, E, and G. In order to complete the ABC classification, Ignatova [16] proposed three voltage dips combined with swells. The H type, in which the main voltage decrement is in two phases and the remaining phase increases its voltage, and the I* and I** type, whose main voltage decrease is in one phase and the remaining phases experience a voltage increase. The equations that represent the three-phase voltage in ABC classification are show in the Appendix.
These ones are obtained by combining the contributions of Bollen [1] and Ignatova [16] and expressing the equations as a function of the same variables used in the SC criterion, i.e., the characteristic voltage V and PN factor F. SC classification is more general than ABC classification, and includes most of the asymmetrical voltage events, which is suitable for a general classification. On the other hand, since ABC classification is more descriptive than CS, it provides more details when a relation between the voltage event and the fault is required. As a consequence, ABC classification is adopted in the three algorithms analyzed in this work as in [17]- [19].

III. Algorithms for Voltage Event Classification
This section presents the operation principles and the main characteristics of the evaluated classification algorithms.

III.1. Symmetrical Components Algorithm (SCA)
This algorithm, proposed in [13], uses the following expression for events classification: where V -1 , V +1 are the negative and positive sequence components of the fundamental voltage, respectively. T is rounded to the nearest integer and can take values between 0 and 5. From this parameter, it is possible to define the kind of dip in the SC classification, as it can be seen in Table I.   TABLE I  DIP TYPE IN THE SC CLASSIFICATION OF SCA (SECOND ROW) AND OF SPA (THIRD ROW) As it was mentioned in the previous section, the ABC classification is more suitable to adequately characterize voltage dips.
Hence, in order to obtain the dip type based on the ABC criterion, it is necessary to obtain three parameters in addition to the T parameter: the characteristic voltage V, the PN factor F and the zero-sequence component of the fundamental voltage V 01 .
The parameters V and F are estimated from (3), presented in [15], and V 01 is obtained from (4):  is a 60° phase shift and V a1 , V b1 and V c1 are the phasors of the three phases.
Finally, using parameter T and values V 01 , F and V in Table II, the dip type in ABC classification and the affected phases are obtained.

III.2. Six-Phase Algorithm (SPA)
This algorithm uses the fundamental RMS single voltages (5) and line-to-line voltages (6), removing the zero-sequence fundamental voltage [13] to classify voltage dips: Voltage dip, according to SC classification [7], is determined by the lowest RMS voltage, as shown in Table I.
Like in the SCA algorithm, it is interesting to know the type of voltage dip based on ABC classification, as well as the affected phases. This is achieved by using the lowest RMS voltage in Table II (third row) and the value V 01 , which is calculated with the expression (4), V and F.
In this case, it can be approximated as follows:

III.3. Space Vector Algorithm (SVA)
In [14], [16]- [17], a technique for detecting, classifying and characterizing voltage events based on the analysis of the space vector trajectory in the complex plane is presented.
The space vector transformation consists in calculating a phasor (8), whose magnitude and phase angle are representative of the three-phase voltages.
Hence, the amount of information required to represent the three-phase voltages is reduced: a phase-shift of 120°. In a symmetrical three-phase system, the three voltages have the same magnitude and are rotated 120º.
Hence, the trajectory of the space vector is a circle of unitary radius.
In [17], it is demonstrated that a symmetrical dip describes a circumference with a radius proportional to the dip depth, while an asymmetrical dip describes an ellipse. To distinguish between different kind of dips and swells, the parameters of the ellipse (9) and the zerosequence voltage (4) are used: where r ma , r me , φ inc and SI are the major axis, minor axis, inclination angle and shape index, respectively, while φ +1 and φ -1 are the phase angle of the positive and negative sequence components of the fundamental voltage, respectively. This classification criterion is provided in Table II, where it can be noticed that, unlike the previous methods, SVA allows to classify not only the voltage dips but also the voltage swells affecting one phase (type H) or two phases (type I * and I ** ).

IV. Voltage Event Classification Issues
The state-of-the-art includes some analysis about the errors in the classification algorithms under certain conditions. However, the causes of this misclassification are not explained in detail.
This work shows that these errors result not only from the particular architecture of each algorithm, but also from the limitations inherent to the classification criteria.

IV.1. Phase Jumps and Asymmetrical Phase Rotations
As described in [25]- [27], the electrical network faults affect the voltage magnitude in the PCC, and produce phase jumps and phase rotations.
Phase jumps can be attributed to two causes.
Firstly, the difference in the X/R ratio of the impedances between PCC and the fault, Zf, and between PCC and the source, Zs, originates phase jumps in the PCC retained voltage.
An exception occurs when Zf >>Zs. Analytical description can be found in Appendix. Secondly, the propagation of unbalanced events through transformers that use asymmetrical connection schemes (for instance, D yn ) can lead to phase jumps in points of the electrical network fed by such transformers, even though the original event does not have such phase jump [26].
Regarding the parameters measured by the classification algorithms, phase jumps produce an argument of V different from zero.
On the other hand, asymmetrical phase rotations originate when the positive and negative sequence impedances between PCC and the source are not equal.
In this case, an asymmetrical phase rotation affects the argument of F. This feature does not occur with passive elements like lines and transformers, but it is possible in the proximity of generators and motors.
Notice that the SC criterion considers that the zero sequence voltage component is zero (V a1 +V b1 +V c1 =0); then, from the three expressions that define the threephase voltage, only two are linearly independent. When this feature is included in the six expressions corresponding to each type of dip (for instance, expression (1) for type Ca), each expression generates a set of two equations with two unknowns, V and F.
Consequently, when analytically solved, this system could make that the same three-phase voltages (V a1 , V b1 , V c1 ) result from any of the six events, each of them with a particular combination of V and F. However, these parameters cannot take any value. In general, |F|≥|V| and the module of V can take any absolute value in so far as |V|<0.9 p.u.. Finally, the phase angle of F or V will be different from zero when phase rotations or jumps occur [15].
In spite of the fact that the amplitude and phase angle of V and F would allow to discard some of the six combinations due to their physical impossibility, the existence of high asymmetrical phase rotations or phase jumps that disturb the three-phase voltages will make it difficult to distinguish one kind of event from the other. This is a major drawback of the classification criterion in general, regardless of the classification algorithm implemented.
This feature is not analyzed in the literature; what is more, misclassification is associated to the classification algorithms performance and not to the classification criterion.
An example of the above mentioned is illustrated in Figs. 2, where the dashed line represents the ideal threephase electric power, the solid black line two ideal events (type Cc (a) and type Da (b)), and the solid blue line the three phase voltages that result from including phase jumps and/or asymmetrical phase rotations in these events.
c, a Note 1: type NS corresponds to the normal state of the network. Notice that an event is defined as ideal when it is obtained only from variations in the three-phase voltages, i.e., there are no phase jumps or phase rotations associated. In the first case, a Cc type event was generated, with V=0.5e -j20º p.u. and F=1p.u. In the second case, a type Da with V=0.54e j10.95 p.u. and F=0.98e -j16.13º p.u. was generated. As it can be observed from this example, there is no way of distinguishing which event occurred, since, in both cases, the disturbed three-phase voltages are the same (blue lines in Figs. 2(a) and (b)), and they were obtained with physically valid combinations of V and F. Therefore, it can be concluded that, regardless of the particular structure of each algorithm, the possibility of a misclassification is always present. Again, this is inherent to the limitations of the classification criterion.
If the same analysis is performed for the ABC criterion, the conclusion is similar, i.e., the same threephase voltages (V a1 , V b1 and V c1 ) can be obtained from different events with particular combinations of phase jumps and asymmetrical phase rotations.
Since this limitation is related to the way in which the classification criteria was proposed, the classification algorithms present the same limitations. A partial study on the errors of SCA and SPA due to phase jumps and asymmetrical phase rotations is given in [13], which concluded that SCA does not yield errors under a phase jump.
The same applies to SPA with asymmetrical phase rotations. These conditions were also evaluated individually, i.e., events with asymmetrical phase rotations and without phase jumps (V=1p.u.), and with phase jumps and no symmetrical phase rotations (F=1p.u.). Since the number of independent variables was limited, the tests did not reveal the drawbacks of the criteria with respect to the event classification. Even though these test conditions allow to analyze the effects separately, they are not representative of an actual operational condition.

IV.2. Limitations of Event Classifications in the Presence of Symmetrical Phase Rotations and Pre-Fault Voltage Variations Regarding Nominal Magnitude
The model used to generate faults with phase jumps and asymmetrical rotations in Appendix is based on the pre-fault voltage and line impedance. In this model, the ratio between positive and negative sequences of the grid voltage is explained by some patterns that are used by the algorithms to classify the event.
A condition that is not considered by the above mentioned model is the symmetrical phase rotation of the grid voltages. This condition can be associated to a switchgear operation, the activation of protection devices that connect/disconnect high power loads, or a fault in the reference used in measurement instruments. For instance, there are situations in which the pre-fault current is high due to the connection of a high inductive load. In these cases, the pre-fault voltage is lower than the source voltage because of a voltage drop in the line, and the phase angle is proportional to the sum of the load and line impedance phase. As a consequence, a symmetrical phase rotation is produced when a fault in the line is located near the load, which may be understood as a reference loss of the three-phase voltages with respect to the prefault voltage. One of the effects of this rotation is an event misclassification by SCA and SVA algorithms, since these algorithms use the module and phase of positive and negative sequence voltages as data for the classification process. However, the classification with the SPA algorithm is not affected since it only uses the amplitude of the three-phase voltages and the three line voltages.
Figs. 3 display a particular example in which an event type Da is synthesized, with a retained voltage V = 0.5e -j20º . Fig. 3(a) shows the event generation under ideal conditions, that is to say, without phase jumps and/or symmetrical rotations, using a pre-fault voltage equal to the unity, and measuring the phase angle of the three voltages respect to the pre-fault phase. Vβ is the phasor whose argument β allows to estimates the value of the parameter T, and with it, SCA estimates the event type.
Under these conditions, the phase angle of the phasor V β = V -1 /(1-V +1 ) , so-called β, is equal to 180°. According to equation (2) and Table I, T=3 and the event is correctly classified by the SCA algorithm. In Fig. 3(b), the same voltage event is affected by a symmetrical phase rotation, which produces a difference of 20º between the three-phase voltages and the corresponding pre-fault voltages. In this case, the positive and negative sequence components are rotated with the same phase angle, but β = 119.01°, and so, T = 2, and considering the criteria of Table I, the event is misclassified as C b rather than as D a .
As it was previously mentioned, the loads may produce variations in the module and phase angle of the pre-fault voltage with respect to the nominal voltage. At the same time, the voltage in each point of the electrical network depends on the load as well as on other factors, such as: -control systems of bulk generators that regulate the voltage level by means of active power control, -the position of voltage regulation tap of transformers, -capacitor banks for power factor correction, -etc.
For the reasons exposed, it is normal that a voltage event in any of the grid phases leads to changes in the non-affected phases. This fact may generate classification errors in the algorithms described in Section III, because these are based on static detection thresholds, which are calculated under the assumption that the voltage in the non-affected phases is equal to the nominal voltage.

V. Performance Evaluation
This section summarizes the results of the performance evaluation of the three algorithms analyzed in this work. The algorithms were implemented in MATLAB®. For the evaluation, the phasors corresponding to each event type were generated based on the equations in Tables A1 and A2 (see Appendix), and the respective fundamental symmetrical components were estimated. The three algorithms were tested with these signals, comparing the outputs to the known input values. Phase jumps and asymmetrical phase rotations were generated by means of variations in V and F module and phase. Symmetrical phase rotations and prefault voltage variations with respect to the nominal magnitude were generated by means of variations in the module and phase of the reference voltage. Note that the first two disturbances are asymmetric, increasing the unbalance of the three phases. While the last two disturbances are symmetrical, affecting all three phases equally.

V.1. Phase Jumps
In order to evaluate the performance of the algorithms under the presence of phase jumps, a set of synthetic dips was generated, according to the ABC classification criterion, with retained voltages from 0 to 1 and phase jumps of the retained voltage ranging from -90° to +90°.  In Fig. 4(a), it can be observed that SCA is not sensitive to phase jumps and classifies all the voltage events type C and D correctly for any combination of retained voltage and phase jump. This can be explained by the fact that the phase between V-1 and (1 -V+1) is not affected by the phase jumps.
Additionally, the simulation of phase jumps only affects the phase of V, maintaining F=1p.u.. The value of F limits the possible combinations of different voltage events, and so the limitations of the classification criteria described in Section II are not evident. However, this test is useful to understand the individual effects. Fig. 4(b) shows the results of SPA algorithm. This algorithm is very sensitive to phase jumps, especially when the dips magnitude is higher than 0.3p.u.. SVA algorithm (Fig. 4(c)) has better performance than SPA algorithm, with a phase jump sensitivity only relevant for dips of magnitudes higher than 0.75p.u. (a higher magnitude of retained voltage involves lower voltage reduction). However, this is not a good characteristic, because dips of reduced magnitude (near 0.9p.u.) are the most frequent in the network.
The results obtained when testing other types of events show some variations in relation to Figs. 4, but the three algorithms have similar performance to that of the classification of dips type C and D.
It is noteworthy that an optimal performance without errors was obtained only for classification of symmetrical events, i.e. A type. This is because symmetrical events have the particular feature that V = F.

V.2. Asymmetrical Phase Rotations
The analysis of the influence of asymmetrical phase rotations on the performance of classification algorithms was conducted by means of two different tests. The first one consisted in the analysis of all the possible combinations of asymmetrical phase rotations and retained voltage for a system in which there are no phase jumps (in other words, the X/R ratio in PCC is equal in the source and the fault side). Even though this seems to be a particular case, the objective is to analyze the performance of the algorithms with regard to the asymmetrical phase rotations, without considering the effects of phase jumps. The results of this analysis for dips type C and D are shown in Figs. 5. The second test consisted in the analysis of all the possible combinations of phase jumps and retained voltages, as previously described in Section V.1, but considering an asymmetrical phase rotation of -20°, in order to understand the effects of the combination of both disturbances. Figs. 6 display the results of this second analysis for dips type C and D.
Despite the fact that Fig. 4(a) shows that SCA algorithm is not sensitive to phase jumps, from Fig. 5(a), it can be noted that this algorithm is very sensitive to asymmetrical phase rotations. As a consequence, when the grid voltage is disturbed by a phase jump and a phase rotation, SCA algorithm fails, as shown in Fig. 6(a).
Basically, this behavior is explained by the variations that asymmetrical rotations introduce in the arg{V -1 /(1-V +1 )}.
The analysis of SPA performance shows that the algorithm is very sensitive to both phase jumps and phase rotations (Figs. 4(b) and 5(b)). This behavior is partially evidenced in [13], and this error is attributed to the unbalance of the three-phase voltages, where the minimum voltage is actually the second minimum voltage.
Finally, the performance of SVA algorithm is similar to that of SCA, that is to say, it is very sensitive to phase rotations (Fig. 5(c)). This behavior is explained by the inclination angle (φ inc ) of the ellipse described by the rotation of the space vector, and giving rise to classification errors. Bands of approximately 30° where SVA algorithm classifies the voltage events correctly can be noticed. These bands are the consequence of the ±15° tolerance bands described in Table II. Finally, Fig. 6(c) illustrates how the algorithm performance degrades notoriously when the phase jumps are combined with a constant asymmetrical phase rotation.

V.3. Symmetrical Phase Rotations
In order to verify the issues described in previous sections, the algorithms were evaluated with symmetrical phase rotations. Fig. 7(a) shows that SCA algorithm fails in almost all cases. On the other hand, SPA algorithm is not affected by symmetrical rotations because they do not change the absolute magnitudes of the six voltages employed ( Fig. 7(b)). Finally, Fig 7(c) lists the results of SVA algorithm, which evidence the classification errors when the rotation angle is greater than the 30° bands derived from the algorithm thresholds.

V.4. Pre-Fault Voltage Deviations
Most references in the bibliography do not consider the case when the pre-fault voltage differs from the nominal or declared voltage (IEEE 1564-2014) and assume ideal operating conditions. To make a more comprehensive study, an evaluation of the classification errors produced by pre-fault voltage deviations with respect to the nominal voltage was carried out. The scenario was the particular case of a symmetrical prefault deviation, in which the three phases were equally affected. Figures 8(a)-(c) show the results of SCA and SPA failure when the pre-fault voltage is lower than 0.967p.u. and 0.96p.u. in the case of SVA. These results were obtained for a balanced system; and errors increase if an asymmetrical condition is considered.

VI. Discussion About the Causes of Algorithm Errors
Below is a brief discussion of the errors in each specific algorithm. The aim is to understand the reasons why the chances of a univocal classification get reduced with these algorithms and criteria.

VI.1. Symmetrical Components Algorithm
Taking into account equation (2), the SCA classification is obtained from the positive and negative sequences components [13].
It fails in the presence of asymmetrical and symmetrical phase rotations because of the formulation of the algorithm itself and the ambiguity of the classification criterion when high phase rotations are considered. In the work by Bollen [13] the algorithm is modified, adding a constant phase angle of 20° to the phase argument estimated with the sequence components (10). This modification mitigates the errors when the phase rotations (symmetrical and asymmetrical) are negative and of little magnitude. This partial solution does not include all the possible cases in actual implementations: The positive and negative sequences should be known in order to understand the reasons why a misclassification occurs when equation (2) is applied.
Then, in this work, the transformation of symmetrical components is applied to the six equations that describe the voltages dips C a (1), C b , C c , D a , D b , and D c , obtaining (11): Vβ is a voltage phasor |V β |e jβ that represents a particular asymmetrical phase rotation for each event type. In other words, the positive sequences as a function of F and V parameters are independent of the event type, while in the case of the negative sequence there is a relation with the event type. From equation (11), it can be deduced that it is possible to characterize each event type with the phasor V β , obtaining two solutions: From (12) and (13), two new equations can be derived to estimate T parameter and, consequently, classify the event without errors in the presence of phase rotations and jumps: The dilemma is that V and F are variables estimated from T V and T F , respectively, so, in order to calculate either of them, it is necessary to implement a recursive or iterative algorithm.

VI.2. Six-Phase Algorithm
This algorithm fails in the presence of big asymmetrical phase rotations and jumps. The reason of this classification error may be ascribed to the minimum RMS voltage that allows to determine the event type, which differs from the voltage prior to the phase jump.
A possible mitigation is to measure the first two minimum RMS voltages, estimate F and V values, and determine which of the possible events is more representative for the values of the characteristic magnitudes.

VI.3. Space Vector Algorithm
The performance analysis of the Space Vector Algorithm indicates that the classification errors occur in the presence of phase jumps associated to voltage dips of reduced magnitude, and phase rotations higher than 15°.
They originate in the deviation of the characteristic angle (φ inc or φ 01 -φ +1 according to the event type) with regard to the event type classification thresholds.
These errors may be partially solved by estimating the phase rotation or jump using parameters V and F.
Nonetheless, as in the previous analysis, the ambiguity of the classification criteria prevents determining the event type without errors.

VI.4. Pre-Fault Voltage Deviation
The classification errors produced by the variations of the pre-fault voltage with regard to the nominal or declared voltage deserve special mention. In this case, it is necessary to work with variables that depend on the pre-fault voltage measurement and not on the nominal magnitudes, allowing an automatic and dynamic adjustment of the classification threshold.
An interesting approach is to estimate the pre-fault voltage in accordance with IEC 61000-4-30 Standard [21], where a sliding reference voltage is used.

VII. Conclusion
This work presents an analysis of the drawbacks and constraints of voltage event classification criteria and algorithms. The analysis was performed based on two criteria and three different algorithms, whose performance was tested under the same operational conditions. The so-called Symmetrical Components and Six-Phase algorithms, originally proposed for SC criteria, were extended to contemplate ABC classification criteria, following suggestions from the literature; and a rigorous performance analysis of SVA, not available in the literature, was introduced.
It was demonstrated that the performance of the three algorithms is limited when voltage events are combined with phase rotations (symmetrical and asymmetrical), phase jumps, and/or pre-fault voltage deviations with regard to the nominal voltage. All these disturbances are common in actual networks and are directly related to the grid and load characteristics.
Despite the fact that some of these errors were reported in the literature for some of the algorithms, this work shows that these errors are not only produced by the algorithm architecture but also by the ambiguity of the classification criteria, which is not treated in the respective literature. In the case of SCA, (14) or (15) are more general and correct equations as compared to (2).
However, they depend on the event voltage type, so a recursive or iterative algorithm should be implemented. SPA and SVA algorithms can also be improved by estimating phase jumps or rotations.
Nonetheless, in the cases where the voltage event is combined with severe phase jumps or rotations, the classification becomes impractical. Based on these results, the need of a new classification criterion arises, using statistical information of actual faults, and more comprehensive description indices.
Finally, the importance of the voltage reference used in algorithms like SCA and SVA is underlined. Taking into account that both algorithms are designed for threephase voltages in relation to the pre-fault voltage, the cases in which the system is not based on this reference should be considered in order to obtain good performance in recursive implementations.
Single-phase faults: where Z f+ and Z S+ are the impedances of the positive sequence grid; Z f− and Z S− are the impedances of the negative sequence grid; while Z f0 and Z S0 are the impedances of the zero sequence grid. Considering the impedance ratio Z f /Z S =λe jα in (1-A), the following expression is obtained: where λ=|Zf/Zs| is a measurement of the electric distance to the fault and α=arg(Zf/Zs) is a measurement of the maximum phase jump that could happen during a voltage event. This arises when analyzing that in expression (5-A), if λ tends to zero, |V * | tends to zero too, and arg(V * ) tends to α. Defining h=λe jα /(λe jα +1), it can be observed that, when (Z f /Z S )>>1, or when the X S /R S ratio is equal to the X f /R f ratio, arg(h)=0 and hence phase jumps do not appear in PCC. Generally, arg(h)≠0 and so phase jumps appear in PCC, whose maximum value could reach α, for instance, when (Z f /Z S )<<1.


extfieldsofstudy:
Physics
Computer Science

text:
Inexact matrix exponential preconditioner for implicitly restarted Arnoldi method in fluid dynamics stability problems for parallel heterogeneous architecture

We wish to solve the problem of linear stability analysis of stationary solutions for a fluid dynamics problem governed by 3D Navier Stokes equations. The discrete problem is of order O(106) - O(107) so we use parallel computations to perform linear analysis. We aim at the heterogeneous computational architecture based on central (CPU) and graphical (GPU) processing units. In the paper we use the implicitly restarted Arnoldi method (IRAM) to recover the most dangerous eigenvalues. The IRAM has poor convergence without specific transformation of the matrix spectrum. Standard methods require either long computations or inverse of the transformed matrix. This can be troublesome in fluid dynamics problems where matrix is not available explicitly and is non-symmetric, ill-conditioned and can be fully filled. We propose the usage of the inexact matrix exponential with shift-inverse transformation using only matrix-vector product. The inverse is performed using Krylov methods using only matrix-vector product. We prove convergence bounds for the GMRES method and show that this type of matrix transformation accelerates GMRES convergence. We demonstrate application of the method in recovering eigenvalues close to the imaginary axis in some model problems on some matrices from the Matrix Market. We also demonstrate the application of the method for some fluid dynamics problems. We demonstrate convergence, acceleration and efficiency on the multiple GPU cluster.

Introduction
During the bifurcation analysis of fluid dynamics problem we need to solve the eigenvalue problem for the most dangerous part of the spectrum at a certain point in parameter space. We assume that by some procedure we derive the finite dimensional system of equations that is solved on a computer: where v(t) ∈ R n is a vector of a discrete solution of the original system, n is large (typically from 1 · 10 6 to 1 · 10 8 ), M is a linear operator (also called mass matrix, can be an identity E) and F is a discrete operator, that is usually combined of linear and nonlinear parts (such operator F is called Ozeen operators for incompressible flow equations), f ∈ R n is a discrete source term and R is the parameter vector. Boundary conditions are included in M and F operators. For 2 some solution v 0 at R 0 of (1) we derive the linearized system for ininitialismal perturbations y = v − v 0 : where A| v 0 is the linear operator of F at solution v 0 and δR is the variation of parameter vector. Operator A can be represented in operator form (as operator-vector application) or in matrix form, in this case A ∈ R n×n . We would like to know the most dangerous part of the spectrum sp(A) that is responsible for the linear stability loss for δR variation. So have two problems of linear stability. The first one is related to periodic and quasi-periodic solutions and in this case the left hand term is not zero. This problem usually poses no computational challenge for eigenvalue solvers. In case of the zero left hand side we interested in the stability of the stationary solution by finding the part of the spectrum closest to the imaginary axis. This problem poses substational computational challenges for the eigenvalue solvers especially in HPC implementation. In order to solve (1) and (2) with high precision i.e. perform Direct Numerical Simulation (DNS), one must use High Performance Computations (HPC). In this papaer we are aiming at the heterogeneous computational architecture that consits of Central Processing Units (CPUs) and Graphical Processing Units (GPUs) in configuration where number of GPUs and CPUs are installed in a single node and these nodes are connected by fast interconnect. This type of HPC cluster configuration is a standard configuration for almost all Top500 machines. We are mostly interested in GPU dominated machines (up to 8 GPUs installed on one node) because GPUs decrease power consumption and machine cost. Let us give a brief survey of methods that are already developed and successfully used during the linear stability analysis of systems (1), where F is the Ozeen discrete operator. The solution of an eigenvalue problem for (1) in case of incompressible flow is conducted using Krylov subspace methods with exponential estimation preconditioner [20] or shift and inverse preconditioner [1]. One of the first such papers was [7] (and some other papers by these authors), where a 3D Floquet stability analysis was applied to a circular cylinder 2D wake stability problem. The fluid flow solver was a spectral element method for 2D equations where Floquet analysis was used through the introduction of infinitesimal three-dimensional perturbations and derivation of a state-transition periodic monodromy matrix. Such approach can be successfully used for the analysis of first bifurcations but the condition number becomes large still for higher Reynolds numbers. Using this approach bifurcation analysis for incompressible flow problems were successfully conducted in papers, for example in [2,3]. Another successful approach is developed by Sanchez and co-authors, some results can be seen in [4,5,6] also using Krylov methods, namely Arnoldi method. The linear stability of a 3D compressible subsonic flow in a cavity is analyzed in [8] using ARPACK implementation of the implicitly restarted Arnoldi method and a sixth-order compact finite-difference scheme for spatial discretization. Other methods exist that can be successfully applied to these problems. For example the novel method based on two sided Newton method [22,23] is successfully appled to find linear stability and optimal disturbances in many fluid dynamics problems. Such approach can be used for HPC application for systems where a basic Newton-Raphson method is tuned and efficiently implemented.
In all cases we can see that there are no methods of eigenvalue analysis that are specifically aimed on multiple GPU architecture. Our goal is to present a method that is used in many our papers to perform eigenvalue analysis on multiple GPU orientated architecture. For the solution of the eigenvalue problem for (2) we use the use Implicitly Restarted Arnoldi Method (IRAM) [20,21,36]. We chose this method for it's efficiency and need only to have a matrix-vector linear operator to execute Arnoldi process. We solve the following eigenvalue problem: where f () is a function that transforms matrix spectrum in some manner in order to accelerate convergence of the IRAM. If IRAM converges, we obtain diagonal matrix f (Λ) and eigenvector matrix V of the desired part of the spectrum (defined separately for particular f ()). Original eigenvalues are recovered, say, by the Galerkin projection: where desired eigenvalues of the original problem are stored in matrix R. The IRAM is based on the construction of the Krylov subspace (K) dimension p: where u 0 is a starting vector with u 0 = 1 and we designate u n := f (A)u n−1 . Every basis vector in (5) is orthogonalized with respect to all other basis vectors by the Arnoldi process. More information on the IRAM is given in [20] and currently used multiple GPU implementation is described in [36]. The method convergence well if desired part of the spectrum is well separated from the rest of the spectrum and is located on the boundary of the convex hall constructed over the spectrum. In other words, by applying f () we wish to transform the desired part of the spectrum to the boundary and separate it. One usually applies the following transformations: Exponential transformation: Polynomial transformation: Shift-inverse transformation: Cayley transformation: Inverse Kronecker sum transformation: Let us assume that we have two leading real eigenvalues λ 1 and λ 2 close to zero with 0 > λ 1 > λ 2 , |λ 2 | < 1. The matrix exponential transformation (6) would be ideal for the IRAM preconditioning only if it was explicitly available. One would seek transformed eigenvalues with maximum magnitude and recover original eigenvalues with maximum real part. However one usually applies time stepping code to calculate (6) where parameter µ = ∆t is a time step in the numerical method. In this case a error on each iteration is multiplied by the factor if e (λ 2 −λ 1 )∆t ∼ 1, where ∆t << 1. Still this method is used to compute eigenvalues in hydrodynamic stability problems, e.g. [2,3], when unstable eigenvalues are sought. A polynomial transformation (7) applies coefficients c j such that a desired part of the spectrum is increased and the rest is damped. This method is used in ARPACK optionally. However it is difficult to apply this transfirmation to separate eigenvalues close to the origin, so, the method is not used in stability problems that we are interested at.
Shift and inverse (8) and Cayley (9) transformations are, probably, the best candidate to precondition eigenvalues located near a specific point (σ = 0 for those near the origin) for the IRAM. The error will be rediced by the small factor λ 1 /λ 2 for |λ 1 | << 1. One, then seeks eigenvalues with maximum magnitude for transformed spectrum and reovers eigenvalues, closest to point σ. Direct inverse in (8) is never used in practice and one has to solve a system of linear equations in the form (A − σE) u n+1 = u n for the IRAM in order to find a new Krylov vector u n+1 in (5). However, it is well known [9,10] that the condition of A can be very poor, especially for high Reynolds numbers and compressible flows. So the inversion in (8) becomes a huge problem that is to be solved in order to solve the eigenvalue problem. The same is true for (10). The condition of such matrix is worse than that for the Jacobi matrix and problem with preconditioners emerges. This problem is worsten for HPC GPU oriented calculations because GPUs have limited amount of device memory and efficient preconditioners for advection-dominated flows (like ILU factorization) require explicit matrix storage. We present the method that allows us to solve the problem of eigenvalue analysis efficiently on HPC multi GPU computational architecture without explicit matrix storage or matrix graph reconstructions.
The papaer is laid out as follows -first we give the preconditioner theory and theoretical convergence estimates. We also include information on the execution algorithm. Then we describe shift-inverse implementation with specific preconditioner for GPU-oriented computional architecture that we use in our calculations. Next section includes testing problems (matrix market examples and fluid dynamics benchmarks) and computational efficiency estimates of the suggested algorithms. Finaly we present some results for 3D fluid dynamics stability problems that were successfully solved by the suggested algorithm. Next, the conclusion follows.
All calculations are performed using an 8-core (16 in hyper-threading) Intel Xeon E5-2640V2 Ivy Bridge and up to six k40 NVIDIA GPUs assembled in one chassis.

Inexact matrix exponent preconditioning for Arnoldi eigensolver
We denote the left half-plane as C −,0 : {∀(a, b) : a + ib, a ≤ 0} and right half-plane as C + : {∀(a, b) : a + ib, a > 0}. Assume that the spectrum sp(A) ⊂ C −,0 . Then we apply the following two stage transformation for the Krylov basis vector (u n+1 ): The first transformation f 1 (A) = (E + µA) m is the inexact exponent transformation (6) which is executed using explicit Euler's method (or another timestepping method used during the solution of (1)) with m timesteps and µ such that these iterations are stable. The contraction is controled by number of iterations m. Eigenvalues λ 1 > λ 2 are transformed to region (1 − , 0) on the complex plane, where = (1 + ∆tλ 1 ) m . If λ 1 = 0, then it is transformed exactly to (1, 0). By this transformation we have sp(f 1 (A)) ⊂ B (1) (0, 0) (unit 2D ball with the center at the origin). The next transformation is a shift in the manner of shift and inverse transformation. The shift σ defines the part of the spectrum to be sought: σ = 1 − + δ. We are interested in the spectrum near the imaginary axis, so δ is selected such that the solution of the linear system (11) convergence rapidly and the spectrum of the resulting mapping is well separated.
To justify our approach may use the following . Let the GMRES method appled to the solution of Bx = b generate residual r q on the q-th step. Then the GMRES method with q steps convergence to the, at least, given precision Proof. Let us consider the q-th residual r q = b − Bx q of the GMRES method [11]: where P q denotes the set of polynomials of degree at most q and with value one at the origin. Following [11] we derive: International Now, by our assumption sp(A) ⊂ C −,0 so mapping (11) will transfer the whole spectrum of the problem to the B 1 (−1 − δ, 0). Analytical complex function reaches it's maximum at the boundary, i.e. max k |p(λ k )| = max Then we can use trial polynomialĈ q and Chebyshev polynomials C q for the boundary S 1 (−1 − δ, 0) [12] and arrive at: Using (14) and (15) we have: followed by (12) after reduction by 2. This completes the proof.
Technically it means that the convergence doesn't explicitly depend on the size of the matrix A. The dependence on the k(V), as it is said in [13] makes the estimate (13) bad. However the bound of k(V) = M is in the logarithm and logarithmically week and dependence on δ is in the base of the logarithm and very intensive. We demonstrate this fact in figure 1.
Since calculation of matrix exponent is a real operation this method can only be applied to pure real eigenvalue detection. To correct this deficiency we introduce the rotation parameter θ = σ[cos(ϕ) + i sin(ϕ)]. Then the transformation (11) is replaced with the following system: where we have additional basis vectors responsible for real u Re and imaginary u Im mapping.
The angle ϕ ∈ [0, ϕ max ], and maximum angle of rotation ϕ max is known due to the following Proposition 2. Let a timestepper that is used to integrate (1) has no dominant stability limit due to diffusion. Let C be a Courant number from the timestepper and let us use the transformation (11). Then: in transformation (17).
Proof. Let us note that µ is selected such, that µC ≤ 1. If the numerical method is stable, then the spectrum of the operator is bounded by the stability region of the method. Maximum is reached on the imaginary part of eigenvalues since diffusion is not a limiting factor. Knowing the C constant we know a bounding region on an imaginary axis. This bounding region is roughly approximated by a bounding box. The most unstable spectrum is located at the neighbourhood of the imaginary axis. Taking all this into account we consider the boundary eigenvalue λ 1 such, that |Re(λ 1 )| << 1, Re(λ 1 ) < 0, |Re(λ 1 )| << |Im(λ 1 )|. Let us apply the transformation (11) . In this case the location of a transformed eigenvalue f (λ 1 ) on the imaginary axis can be bounded as: So the maximum angle of rotation ϕ max that brings this corner eigenvalue to the real axis by transformation (17) is bounded by (19) from bellow and by π from above.  Figure 1. Illustration of the q vs δ dependence by (12) for a matrix with M = 1 · 10 12 and ε = 1 · 10 11 .
The system increases its size twofold resulting in additional storage space which is a significant drawback. Another drawback results in the introduction of rotational grid that must be used in order to find all leading eigenvalues with nonzero imaginary parts. However this approach, in authors experience, is more stable in practice than the Cayley transformation. The convergence properties don't change since the mapping (17) has no affect on k(V).
In order to recover the original eigenvalues of the problem from the mapped operator (17) we perform the following projections: Desired eigenvalues of the original problem are contained in the matrix R. Notice, that pure real eigenvalues are doubled by this projection.
Let now there exist the part of the spectrum sp 1 (A) : sp(A) ⊃ sp 1 (A) ⊆ C + with leading eigenvalues Re(λ 1 ) > Re(λ 2 ) > 0. It means that the main solution becomes unstable. Then we apply inexact exponential transformation (6) and eigenvalues are transformed as (1 + ∆tλ j ) m , j = {1, 2}. So we controll the separation of the desired spectrum part sp 1 (A) by increasing ∆t or increasing m. This approach for multiple GPU architecture and stability analysis is illustrateded in [36]. For the linear system solution we use restarted GMRES(n) or BiCGstab(L) methods (the latter has better convergence in practice). Schematically, the described algorithm for inexact matrix exponent preconditioner is presented in figure 2.
From this figure we can estimate number of matrix-vector multiplications as: where K IRAM is the number of vectors in Krylov subspace of the IRAM, I Arnoldi is the number of Arnoldi iterations, K GM RES is the number of basis vectors in GMRES, I GM RES is the number of restarts for GMRES method and m is number of iterations in the inexact exponent mapping (11). Typical value of op(M V ) is about 60000-80000 from our experience. This is compareable with the typical M V count for the linear system solution during the standard shift and inverse process with good preconditioning. However, this approach is way more efficient on the multiple  Figure 2. Schematic discription of the inexact exponent shift and inverse IRAM algorithm GPU systems since minimal PCIExpress transfer is used during iterations and only matrix vector applications are used by calling linearised timestepper for the given vector.

Shift-inverse preconditioning on multiple GPUs
We also implemented shift-inverse and Cayley transformations (8) for multiple GPU architecture using simplified approximate preconditioner for the solution of the linear system.
Survey and analysis of preconditioners for fluid dynamics problem or GPU architecture is given in [14,15,16] and for shift inverse methods -in [24,25]. One can observe that the iterative convergence properties of the matrix f (A) are worsten by the decrease of diagonal elements in (8). The preconditioned shift-inverse method is solved iteratively using Krylov-type method (GMRES or BiCGStab(L)). Let x j be an approximation of the eigenvector and B := (A − σE). An iteration of the Krylov method can be formulated with the right preconditioner P as: Olshanskii2012 We construct the preconditioner by using the linearisation on more simple operator with close spectral properties (pseudo-compressibility method for incompressible flows or first order method for compressible flows) with simple Lax-Friedrichs fluxes and simple diffusion discrete operator (7 point Laplace operator in 3D). It is easy to show that the these methods generate diagonally dominant matrices, even for large Reynolds numbers. The Jacobi matrix is found using analytical differentiation of the Lax-Friedrichs flux. Then the matrix can be inverseed by the application of the Block Jacobi or Block Symmetric Gauss-Seidel methods and applied as the preconditioner that gives a rough estimate of the linear operator spectrum. We use Geometric Muligrid (GMG) method [17,18,19,32] to accelerate convergence of the preconditioner. It was tested in many applications on GPU architecture and shown good parallel efficiency.
Let us denote GM G (A, b) as a GMG set of V or W cycles for the solution of the linear system with matrix A and the RHS b. Now we construct the preconditioner. First we split our operator F(u) = N(u, u) + Cu to nonlinear part N() and linear part C. The nonlinear part is found using Lax-Friedrichs flux. Then we apply linearisation that gives us linear operator P| u 0 := N u 0 + C near the solution u 0 . The linearisation of N() and linear matrix C becomes split in diagonal and off-diagonal parts P := L + D + U: Here λ max is the spectral radius of the Jacobi matrix that can be estimated by the maximum absolute value of the spectrum in the linearized system. Observe that matrices L, U and D consist  (1) dimension. For gas dynamics problems these are blocks of 5 × 5 and 3 × 3 for incompressible flows. First we add shift to the system by the scalar parameter σ that influences only diagonal block matrix D, obtaining modified diagonalD = D − σE.
Then for each smoother in the GM G(P, r) we perform the following sweeps: In this case Krylov residuals r n+1 are constructed using preconditioner: Sweeps (24) are performed in the menner of LU-SGS method analogous to the suggested approach in [26,27]. First we decompose our computational domain into separate subdomains for each computational device (GPU or CPU) and color these domains in such manner, that adjusted boundary elements belong to domains with different coloring, say "blue" and "green". Next, we select internal and boundary elements for these domains. First, we perform sweeps in "blue" internal cells, then we perform sweeps in half of "green" internal and all boundary cells, then we perform sweeps in "blue" boundary cells and, finaly, use the rest of the "green" cells. Inside every cell subdomain we use red-black coloring to perform GPU or OpenMP parallelization. Boundary cells are exchanged using CUDA-aware MPI. After that we perform either restriction or prolongation on GMG and perform sweeps again on different grids. It is noted in [24] that a standard preconditioner to B is not very efficient so one must apply a tuned preconditioner formed by Sherman-Morrison formula. However we were unable to find any implementation of such preconditioner for HPC architecture and this problem is the topic of another research.

Testing problems and efficiency
In this section we give results on testing problems, efficiency estimates and examples of the method application.

Matrix Market examples
First we present some test on standard matrices. Small matrix rdb450 from matrix market collection [31] is used to illustrate the application of the method. This non-symmetric matrix is originated from the discretization of reaction-diffusion brusselator model. Different stages of the method application are illustrated in figure 3. Green and red circles represent eigenvalues that are found by the application of the described inexact matrix exponent preconditioner and IRAM. We seek least stable part of the spectrum, i.e. eigenvalues with maximum real parts. Original spectrum of the matrix and desired eigenvalues are depicted in figure 3a. First we perform inexact matrix exponent transformation, result is presented in figure 3b. Observe, that the desired eigenvalues are now largest magnitude eigenvalues (LM) and are located close to the line Re(z) = 1. Then we perform shift with σ = 1.05 transforming the whole spectrum close to the line Re(z) = 0.05. Then we inverse the transformed matrix using GMRES method (in this case it converges to machine epsilon with 4 iterations) and the rest of the spectrum is inverseed analytically for illustrative purpose. The result is presented in figure 3c. We can observe that the method successfully recovered eigenvalues near the origin. However, one more complex-conjugate eigenvalue is not recovered (marked with red circles in figure 3c). In order to recover this eigenvalue we apply rotation to the matrix transformation using (17) with rotation angle 0 < ϕ < π/8 and obtain missed eigenvalues at ϕ = π/11, shown in figure 3d with the transformed spectrum. Next we consider the stress test problem of finding four eigenvalues with largest real part for non-symmetric matrix S80PI n1, sized 4025 × 4025. This matrix has two complex-conjugate desired eigenvalues λ 1,2 = 2.0 · 10 −12 ± 1.87i and λ 3,4 = 2.0 · 10 −12 ± 1.68i. We compare the suggested transformation (11) with Cayley transformation (9). For this test we chose σ = 0.01, σ 0 = −0.01, 0 < ϕ < π/3 with step dϕ = π/50 and use m = 3 with µ = 0.3. We add 10 additional Krylov vectors for filtering in the IRAM, thus Krylov subspace dimension is 14.
Results for Cayley transformation are provided in figure 4. One can see that it is impossible to recover desired spectrum with this transformation, despite good convergence of Ritz eigenvalues. The maped spectrum (middle figure) is tuning Ritz eigenvectors to the boundary. On the other hand we are able to recover desired spectrum using suggested transformation, see figure 5. One can see, that the maped spectrum of the problem is transformed in such way, that the IRAM recovers selected eigenvalues very efficiently. Red circles conrespond to ϕ = π/6.5 and cyan circle -to ϕ = π/7.5. This process takes some time, because we need to solve the eigenvalue problem for each value fo ϕ in accordance with the introduced step size dϕ.
Provided examples show that the suggested preconditioning method can be used to recover the most dangerous part of the spectrum from the linear stability point of view.

Fluid dynamics examples
Now let us consider the acceleration achieved for multiple GPU architecture and compare acceleration and efficiency with shift and inverse transformation using preconditioners (25).
Results of convergence for different methods are provided in table 1. One can see, that the choice of parameters σ, m and µ can heavily influence the method. For example for the same problem size (say 6.4E6) we can obtain convergence for the restarted GMRES method from 1100 to 3600 iterations. unfortunately, there's no a priori recipe for selecting these parameters. This we leave for future work. Second, we can see, that the BiCGStav(L) method converges more rapidly, so in practice we use this method more often.
Verification of the method on fluid dynamics applications is performed in all problems we solved, e.g. [33], [34], [35]. For example we present here the comparison of growth rates as   [35] are presented in figure 6 and wall time for compressible flow problems [34] are presented in figure 7. We can see, that only four iterations of the full IRAM are required to achieve desired tolerance. The same rapid convergence was observed in figure 5. This is achieved by adjusting regulating parameters µ and σ for the particular problem.  We noticed that the computational cost of the suggested method is O(N ) for compressible flow problem and O(N log 2 N ) for incompressible problem using spectral method. It means that the suggested method doesn't introduce additional difficulties during computation compared to the original time stepping method proportional to the problem size. Additionally one can see that the initial stagnation is related to the incomplete GPU loading.
Comparison of the preconditioned shift and inverse method with the suggested method is shown in right figure 7 and done on the compressible problem of the Kelvin-Helmholtz instability [34]. We can notice that preconditioned shift and inverse method is more efficient for small matrices, however, starting with problem size around 1 · 10 6 the presented method becomes more efficient. This can partially be explained by less efficient preconditioner for the shift  Figure 6. Convergence of full IRAM iterations and wall time required to solve the problem of finding 6 largest real eigenvalues in spectral method for incompressible flow for the residual target 1 · 10 −10 in L ∞ -norm. Figure 7. Wall time required to solve the problem of finding 6 largest real eigenvalues using IRAM for the residual target 1 · 10 −10 in L ∞ -norm and ratio of wall time for shift -inverse method using Lax-Friedrichs and GMG preconditioner. and inverse method. One can also check that the efficiency of parallel implementation for the preconditioned shift and inverse method becomes lower with the increase of parallel nodes. This is due to the increase of PCIExpress communications.
Convergence of eigenvalues to true ones in the physical problem is the function of problem discretization. In figure 8 we show the convergence of first four leading eigenvalues to the true ones in the Kelvin-Helmholtz problem as function of domain discretization and show modulo of the first unstable eigenvector.

Conclusion
In the paper we present the method of preconditioning for the linear stability eigenvalue problem. The method is aimed at the implicitly restarted Arnoldi method and can be easily added to the existing linear analysis code for fluid dynamics problems. The main purpose of the method is to be used for stability analysis of stationary solutions on cluster based on multiple GPU parallel architecture. We show that the method has analytical convergence bound and application examples confirm this bound. We also demonstrated efficient parallel application of the method in fluid dynamics problems sized up to 2 · 10 8 elements. Our future work is to tune the method and give suggestions for the selection of tuning parameters. The method implementation is available freely under GPL licence by request to the author.