Transfer Function with Nonlinear Characteristics Deﬁnition Based on Multidimensional Laplace Transform and its Application to Forced Response Power Systems

: In this research, the concept of nonlinear transfer function with nonlinear characteristics is introduced through the multidimensional Laplace transform and modal series (MS) method. The method of modal series is applied to the power systems dynamics analysis in order to consider nonlinear oscillations and modal interactions, which contribute to the response of the system’s dynamic following disturbances. The method of MS allows the inclusion of input excitation functions obtained as Laplace domain kernels superposed to obtain a transfer function. Applying the Volterra series expansion through kernels decomposition, a transfer function with nonlinear characteristics is obtained which incorporates some of the main modal characteristics of the nonlinear system. Following the same schematic procedure, it is possible to determine second and higher order transfer functions. Once the transfer functions both linear and with nonlinear characteristics are determined, a time domain and frequency response analyses can be performed. The methodology is exempliﬁed by denoting the numerical and analytical properties with the application to a synchronous machine-inﬁnite busbar test power system and to a three synchronous machines–nine buses test power system. Bode and Nyquist analysis are utilized to demonstrate the transfer functions accuracy and frequency response.


Introduction
The transfer function is a concept normally utilized in control systems defined as the relationship between the Laplace transform of an output variable with respect to the Laplace transform of an input variable [1,2]. On the other side, practical applications of higher order models in different areas of practice engineering have been researched covering novel specialties [3], i.e., integrated circuits and micro-systems, civil engineering, aero-space, vibration problems, microelectromechanical systems (MEMS), among some others. Of course, in the field of practical engineering applications of nonlinear and higher order models are the power systems and some of the problems presented in their operation, such as nonlinear oscillations. The concept of transfer function has been extended to the study of power systems with many contributions to related issues such as matching nonlinear where A is the Jacobian matrix (first order partial derivatives) Since Equation (2) is a linearized but coupled system, a linear transformation like Jordan canonical form is used in order to express a new system with uncoupled linear and nonlinear characteristics. If x = Uy where U are the right eigenvectors from the state matrix A in Equation (2), the Jordan canonical form system is, . y j = λ j y j + The notation of u P p is the element at pth row and Pth column of U. The last term of Equation (3) represents the control variables. This term is reduced tob ji u i when a single input-single output (SISO) system is under study; on the other side, the term is the summation assuming a multiple input-multiple output (MIMO) system of number m.

Time Response under Impulse Function Input
As it was explained above, the method can be generalized to incorporate the effect of arbitrary excitations on the dynamic system. Since the output function is the result of each individual term from the Taylor series and Jordan canonical form, the complete function has the form in accordance with a Volterra series, that is [27], . y j = y 1 j + y 2 j + y 3 j + . . . + y n j +b j u j .
The system now is expressed by considering the individual and uncoupled terms. Please observe that for this case, the input function has influence in the first order term due to the nature of the exciting function. However, the second, third and higher order are linked with the input function resulting in an influence over the rest of the terms. Additionally, this result explains the modal interaction between linear and higher order terms. Thus, . y j = λ j y j +b j u j , . y 2 j = λ j y 2 j + n k=1 n l=1 C j kl y 1 k y 1 l , . y 3 j = λ j y 3 j + n k=1 n l=1 C j kl y 2 k y l + y k y 2 l + n p=1 n q=1 n r=1 D j pqr y 1 p y 1 q y 1 r . (9) . . .
Applying the Laplace transform to Equations (7)-(9), we have, where To illustrate the proposed method, consider the case of an impulse function u j (t), u j (t) = δ(t) ⇒ U(s) = 1.
From Equation (10), Obtaining the time domain solution by applying an inverse Laplace transform results in, Following the same procedure, it is possible to obtain the Laplace domain kernel of second order term, i.e., When multidimensional Laplace domain kernels N 2 (s 1 , s 2 ) are obtained, it is possible to calculate a single Laplace domain kernel applying the procedure known as association of variables [29,31], that is, . (20) Calculating the respective coefficients, the second order kernel expressed in single Laplace domain has the form, which represents the time domain solution given by, Equations (15) and (19) can be utilized for both cases of time domain analysis and for finding higher order transfer functions. Finally, the complete time domain solution of Equation (6) which is defined in terms of canonical form variables is expressed by, and converted to the original state variables,

Synthetic Example
To illustrate the application of the method, let us assume a synthetic example taken from Reference [32]. A first order nonlinear differential equation of the form, is solved when the input function u(t) is an impulse function, as a first step, and afterwards, a step function is also studied; both using the modal series method. The differential Equation (25) has first and third order elements, which are considered in the final solution. Thus, assuming an impulse function, the equation takes the form, the first order term is calculated as, . Applying the inverse Laplace transform, This particular system does not have second order term. Hence, the third order term is expressed as, where, Applying theorems of association of variables [9], the solution in terms of a single Laplace domain is, Finally, the full solution to the differential Equation (25) for an input impulse function has the form, The time domain solution expressed by Equation (30) is like that described in Reference [32] under the same conditions (input impulse function). This sample has demonstrated the viability of the extension of the modal series method to the inclusion of a forced input signal. Similarly to the previous case, a step input function was analyzed for the same nonlinear differential equation. This differential equation takes a similar form to Equation (25); however, the first order terms, which include the step function have the form, with y(0) = 0. Therefore, Executing the inverse Laplace transform, The third order terms are those of Equation (28), but for the step function takes the form, This equation may be solved using the theorems of association of variables on the multidimensional Laplace transform, described in Reference [29]. Applying the real convolution theorem [31], and Finally, the third order terms in the time domain solution have the form, with the full solution, Equation (38) agrees again with that obtained in Reference [32] for the case of a step input signal. The analytical example has demonstrated the application of modal series methods based on multidimensional Laplace to the differential Equation (25) and their comparison against the solution obtained with a different technique such as the Laplace-Borel transform [21].

Introduction
The transfer function concept has been deeply analyzed as a powerful tool for either invariant linear time or time varying systems, with solid applications to control systems. The concept of the transfer function itself addresses the behavior between inputs and outputs of the dynamic system defined in terms of the Laplace domain. The transfer function concepts to power systems have diverse potential applications, since there are many actions of control which a real power system has to implement and perform. Some characteristics of transfer functions based on power system may be observed [13]: The order of any transfer function is as high as the number of state variables included in the model; i.e., in large power systems there are many thousands of states; most of those poles or eigenvalues are not observable from any signal in the system. An input signal applied to a power system will usually excite many modes, which will be reflected in different output signals and locations in the system. Despite a transfer function is usually defined as a ratio of polynomials, it can be expressed as a sum of residues over a first order pole. Transfer function analysis assumes that initial conditions are out of concern and the output is only a function of the forced function.
The concept of the computation of dominant poles of power system transfer functions has been addressed by some authors [5][6][7][8][9][10]. The concept arises over the consideration of a transfer function in which the poles and zeros are known, determining the main modes that are taking part in the frequency response to the oscillations in the power system. An algorithm based on Rayleigh iterations has been proposed, which exploits its characteristics of numerical properties of global cubic convergence [5]. The participation factors and the applications of such an algorithm are marked when a large-scale power system is studied. An earlier work proposed the use of Prony analysis to determine the residues and with them, to obtain the transfer functions for the design of power systems stabilizers (PSS) applied to a multimachine power system [33]. In this work, the usefulness of transfer function concept arose from determining the frequencies during oscillations in the power system that eventually may have been damped by PSS devices or another flexible AC transmission systems (FACTS) component.

Theoretical Basis of Nonlinear Characteristic in Transfer Function
Halás states in his contribution [19]: "everybody knows that a nonlinear system has no transfer functions". This affirmation is completely valid if we recall that the concept of the transfer function itself is only related to linear systems. However, many techniques like those mainly described in Section 1, arise with nonlinear systems that after some mathematical methodologies allow us to obtain the so-called nonlinear transfer functions. The extension of modal series proposed in this contribution took advantage of this reasoning.
Transfer functions of nonlinear systems must satisfy many properties that are expected from transfer functions [3]. Some of the main properties were resumed here:

•
They characterize a nonlinear system uniquely. Each nonlinear system has a unique transfer function, no matter what state space realization one starts with.

•
They provide an input-output description of the nonlinear system. • They allow the use of transfer function algebra, to combine systems in series, parallel or feedback connection.
The pseudo-linear algebra concept that unified the formalism of algebra related with both linear and discrete nonlinear systems is very important [34].
Assuming a transfer function defined for a linear system, we have, The time response of Y(s) to an impulse disturbance applied in the input function U is equal to the inverse Laplace transform of H(s) considering zero initial conditions in all system states [23]. Thus, where λ i are the poles of H(s) and R i the associated residues. The time response of H(s) can be directly obtained by integrating in the time domain the impulse response of H(s), that is, This case is very illustrative since the time response contains two parts due to the generic disturbance [31]:

•
The first part is called forced response, which consists of the steady state part.

•
The second one, called natural response, consists of the transient part, formed by a sum of exponentials, whose values depend on the applied disturbance.
The natural characteristic response depends on the location of the system poles in the complex plane [23]. Therefore, it is important to point out that modal analysis consists on the computation of these poles and on the determination of their nature and sensitivities [33].
For instance, the transfer function of a single input-single output (SISO) system has the form [7], whose transfer function is, where the eigenvalues are the poles of A and the corresponding right and left eigenvectors are given by λ j , x j , v j . The transfer function can be expressed as a sum of residues R j over first order poles, that is, where the residues R j are R j = x T j C v * j B . A pole λ j that corresponds to a residue R j with large R j /Re λ j is called a dominant pole, which is observable and controllable in the transfer function [8]. Several efforts oriented to research the dominant poles in large power systems have been focused on the transfer function concepts. Almost every developed approach is based on numerical algorithms that allow the determination of poles that have the biggest participation in power system oscillations [5][6][7][8][9][10]. Their extension to frequency analysis is straightforward, since functional kernels expressed in the Laplace domain are easily shifted to the frequency domain. Both SISO and MIMO systems are described by each contribution following the properties of linear transfer functions. The calculation of dominant poles was not be addressed in this paper. So far, only the theoretical basis of the nonlinear transfer function was of concern, leaving room for future developments in this research topic.

Volterra Functional Expansion
The input-output behavior of a dynamic system of the form, may be represented by means of a generalized convolution integrals series [35]. A general convolution integral of order k is defined as follows: let (i k , . . . , i 1 ) be a multi-index of length k with i k , . . . , i 1 of the set {1, . . . , m} and if u 1 , . . . , u m are real value continuous functions defined on [0, T]; the generalized convolution integral with kernel ω i k , ... ,i 1 is defined as [30], for 0 ≤ t ≤ T. After some manipulations, it can be demonstrated that the series written as, is absolutely and uniformly convergent, and is called a Volterra series expansion [35,36]. The original application of the Volterra function to the analysis of nonlinear circuits is due to Wiener [36,37]. Wiener established that an output function y(t) of a nonlinear system is function of its input u(t) and that the two functions can be related by a functional series. The first few terms (up to order three) of the functional expansion are, The nth order kernel of Equation (48) h n (τ 1 , τ 2 , . . . , τ n ) can be called a nonlinear impulse response of order n [36]. Its Fourier transform can be called the nonlinear transfer function of order n. Hence considering the input-output relation it can be written in the form A system H of the form can be broken into a parallel combination of systems H i as shown in Figure 1. The interpretation of such a figure can be resumed as [38]: • The generalized functional representation shows a nonlinear system as a parallel bank of systems H i that are nth order nonlinear systems.

•
They have an impulse-response function h n (t 1 , . . . , t n ) associated with them.
The nth order kernel of Equation (48) ( ) 1 2 , , , n n h τ τ τ  can be called a nonlinear impulse response of order n [36]. Its Fourier transform can be called the nonlinear transfer function of order n. Hence considering the input-output relation it can be written in the form  Figure 1. The interpretation of such a figure can be resumed as [38]: • The generalized functional representation shows a nonlinear system as a parallel bank of systems i H that are nth order nonlinear systems.
• They have an impulse-response function associated with them. Following the same reasoning, a nonlinear system of the form, can be split-up into a sequence of components connected in parallel by the method of Volterra functional expansion [36]. The process is illustrated by Figure 2. For this case, the first component is linear [39], that is, Here, the procedure based on the multidimensional Laplace transform can be used. In this way the linear case can be transformed as, Following the same reasoning, a nonlinear system of the form, can be split-up into a sequence of components connected in parallel by the method of Volterra functional expansion [36]. The process is illustrated by Figure 2. For this case, the first component is linear [39], that is,  Here, the procedure based on the multidimensional Laplace transform can be used. In this way the linear case can be transformed as, The second component is of quadratic nature, i.e., In order to use the theory of the multidimensional Laplace transform, it is necessary to artificially introduce t 1 and t 2 , i.e., and then, Formally, at least F 2 (s 1 , s 2 ) can be inverted to obtain f 2 (t 1 , t 2 ) for which f 2 (t) is the desired output [38]. The method is generalized to higher order cases, hence the third is a cubic component; its output is, which in the Laplace transform domain is, The procedure of taking a number of variables t 1 , t 2 , . . . , t n as equal is called association of variables. The great significance of making the associations in the transform domain lies on the fact that these associations can be made by inspection; this application being a class of problem adapted to this theory. Hence, the total system output is given by, assuming a finite number of terms represented as an nth order nonlinear approximation. It must be pointed out that terms associated to Volterra representation are denominated impulse-response Volterra kernels [31]. Such kernels can be compared with those obtained by the modal series. Following the same reasoning, it is possible to represent a forced nonlinear system decomposed by the modal series method, in the same way as the Volterra framework. By comparing each kernel, similarities from a qualitative point of view for both kernels' sets can be observed. Also, it is important to mention that one of the main applications of multidimensional Laplace transforms centers on the closed form solution of Volterra kernels.

Analytical Deduction
Based on the modal series approach developed above, and considering the definition of the transfer function, which remarks the absence of initial conditions on the transfer function concept, it yields, Following the procedure of association of variables which expresses the referred second and third order terms in a single Laplace variable, we obtain, It can be observed that the poles of the first order transfer function are given by the eigenvalues obtained from the state matrix through Taylor series expansion of the nonlinear system. However, second and third order term poles, besides the poles due to eigenvalues, are due to the combination of two and three eigenvalues, respectively. Figure 3 shows in schematic form the linear and nonlinear flow chart to obtain the closed form transfer function based on the modal series approach.

Classical Model
For the sake of exemplification, a simple power system consisting on a synchronous machine-infinite busbar power system was analyzed. The system was based on classical modeling [40], with extended representation up to third order and assuming a single input-single output. The system has the form, which includes only one transfer function given by, Applying the modal series, it yields, • First order terms • Second order terms which is associated to, • Third order terms This process can be applied repeatedly up to the terms suitable to be considered. In this case, terms up to third order were solved, hence the transfer function with nonlinear characteristics was represented by linear, second and third nonlinear terms. Figure 4 shows the Bode graphics for the transfer functions, and corresponding to the first, second and third order, respectively, associated to the test power system considered in this experiment. Bode analysis was performed by using MATLAB ® control toolbox (R2016a, The MathWorks, Inc, Natick, MA, USA), assuming individual transfer functions. The frequency contributions of second and third order terms in the gain of transfer function were observed.  Agreeing to the reasoning of residues, the corresponding poles (modal combination) referred to 5 the resultant residues are shown in Figure 13. This analysis could be compared with the proposal of 6 dominant poles calculation [20,25] where the main residues and their corresponding poles are 7 obtained using an iterative process. With respect to the complete poles of the transfer functions, the 8 main frequencies detected in the wide range spectrum were 4. 16  The Bode graphs of Figure 4 show the frequency contributions of both linear and nonlinear transfer functions. Of course, the linear frequency transfer function was due to the poles of the oscillatory eigenvalues with frequency f = 1.1096 Hz, such as it is described in Table 1. Furthermore, the frequencies due to the second order transfer functions were caused by the combination of oscillatory modes, resulting in two main frequencies contributions, also detailed in Table 1.  This Table specifies the residue, pole, frequency and damping ratio associated to each mode for the case of the first order transfer function, and the resulting modal combination for the case of the second order transfer function. A frequency of 2.2191 Hz (which can be observed in the Bode graph of Figure 4) was the result of modal combination. Also, to be noted was the substantial increase on residual values of the second order transfer function, dependent of nonlinear coefficients.
Another schematic form to express the frequency response of the linear and nonlinear transfer functions was the application of the Nyquist plots or polar plots shown in Figure 5. In this case study, the mechanical power input, as mentioned before, was applied as the input signal to the transfer function. Figure 5a shows the polar plot for the linear transfer function. According to the theory of Nyquist diagrams, the frequency evolution demonstrated that the system was stable for a defined perturbation, since the polar diagram is not circling the axis in one and the poles are in the left half plane (negative real part). The analysis was made considering transfer functions in open loop.    The waveforms of Figure 5 show the time evolution of each transfer function with respect to their individual residues. Of note are the values of second order transfer function residues that were not considered when only the linear approximation is represented.
Also, the same reasoning was used for the second order transfer function. A zoomed vision of the diagram of Figure 5 is shown in Figure 6, where the frequency of 2.2 Hz, which was contributed by second order modal terms, is observed. In the same way as the linear case, referring to Figure 6, the trajectory obtained by the Nyquist plot is apparently circling at point −1, but in counter-clock direction; this is one condition to demonstrate that the system was stable for the transfer function performed under the parameter conditions given. Nyquist theorems must to be explored and extended in future contributions where control dynamics is incorporated. The waveforms of Figure 5 show the time evolution of each transfer function with respect to their individual residues. Of note are the values of second order transfer function residues that were not considered when only the linear approximation is represented.
Also, the same reasoning was used for the second order transfer function. A zoomed vision of the diagram of Figure 5 is shown in Figure 6, where the frequency of 2.2 Hz, which was contributed by second order modal terms, is observed. In the same way as the linear case, referring to Figures 6, the trajectory obtained by the Nyquist plot is apparently circling at point -1, but in counter-clock direction; this is one condition to demonstrate that the system was stable for the transfer function performed under the parameter conditions given. Nyquist theorems must to be explored and extended in future contributions where control dynamics is incorporated.

Step Input Response
The transfer functions previously defined are subjected to the application of a unit step function of the form [2],

Step Input Response
The transfer functions previously defined are subjected to the application of a unit step function of the form [2], f (t) = 0 for t < 0 f (t) = 1 for t > 0 The unit step function is undefined at t = 0. In the Laplace domain it is defined as, Taking the transfer functions Equations (63) and (64) defined for linear and second order, respectively, and applying the input unit step function, the closed form expressions in the Laplace domain are, ∆ ω 2 (s) =ĉ 1 Now, applying the inverse Laplace transform, the closed form solutions in time domain are obtained as follows, Performing simulations through the time domain expressions of Equations (54) and (55) and according to the assumed parameter conditions for this case study, the waveforms of Figure 7 are obtained for the case of input step function response.
Low frequency oscillations are observed in Figure 7, where the second order evolution was considered. The graphs also show that the system remained stable under this input signal, thus oscillating to eventually disappear. An additional exercise could be performed in this case study, which assumed various damping coefficients, resulting in different oscillation conditions. Figure 8a It was of note that the oscillations varied only in amplitude, keeping the frequency unchanged; as it was expected. For higher values of damping conditions, oscillations disappeared more rapidly than low damping constants.

Three Synchronous Machine-Nine Buses Test Power System
The application of the proposed method was illustrated on a three machines, nine bus test power system. The generation and network parameters, load data and the system operating conditions were taken from Reference [41] and given in Appendix A. The power system was represented by a fourth-order model including the automatic voltage regulator (AVR) representation.

Nonlinear Transfer Functions
Taking into account the modeling of a multimachine power system, and considering the fourth order power network representation, the system was linearized by considering the state variables of each generator and the input variables associated to the power system as the mechanical power inputs and the reference voltages of the automatic voltage regulator, i.e., Hence, the complete model including the power system under study has the form, where F 1 (x), F 2 (x) and F 3 (x) are the Jacobian, Hessian and third order terms matrices [27], respectively. The input matrix B has the form, which for the assumed power system is a (15 × 6) order matrix. For this case study, the transfer functions could be described by choosing the rotor speed deviations as the set of output variables, which resulted in a multivariable transfer function with six inputs-three outputs, in the form shown by Figure 9. Since the system under study has multiple inputs-multiple outputs (MIMO system) it is necessary to incorporate the concept of superposition. Although this is a concept only valid for linear systems, the linearity of inputs and the nature of each kernel represented in terms of the multidimensional Laplace domain could be assumed. For the sake of simplicity, each transfer function was calculated as a single input-single output system thus considering only the effect on the output variable due to a single input variable.
Applying the modal series method, detailed transfer functions including second and third order terms can be described. Following the modal series reasoning, the complete system can be modeled according to linear rules, which in its own can be expressed according to Jordan canonical form, i.e., which for the assumed power system is a (15 × 6) order matrix.
• First order transfer functions: The transfer function matrix has the form, where, s − λ 15 , . . .
• Second order transfer functions: To formulate the second order transfer function, it is necessary to consider the system as a single input-single output in order to form the transfer function matrix, that is, where . . . , n; n = 15; q = 1, 2, . . . , r; r = 6. From Equation (79), considering the selected output variables, which for this case study were the rotor speed deviations (ω 1 , ω 2 , ω 3 ), the individual transfer functions are obtained and grouped to form the transfer function matrix. For instance, the first transfer function is formed as, Energies 2019, 12, 4061 23 of 32 In the same way, the rest of the transfer functions are defined and then grouped in matrix form in order to complete the transfer function matrix as, where , . . .
In the general case, with, i = 1, . . . , o, o = 3; j = 1, . . . , r, r = 6. The input signal is chosen as the input mechanical torque, for which it is possible to determine the transfer function of any state variable. Upon rearranging terms, we obtain the transfer functions. These equations are amenable to multidimensional Laplace transform analysis. We emphasize that, unlike previous approaches, the above formulations allow the study of the system behavior in both the time and frequency domains.
Finally, the complete transfer function for this MIMO system is expressed in matrix form as, • Numerical Analysis The example was conducted following the theory of nonlinear transfer functions through the modal series described above. In this study case, the linear transfer functions were defined by Equation (64)   ❖ Nonlinear transfer functions Now, the nonlinear transfer functions are analyzed. Of note is Figure 11, which denotes the Bode and Nyquist graphics for the transfer function Ω 2 (1,2) (s). It is important to notice the several frequency contributions that the Bode graph shows. This transfer function is an element of the nonlinear transfer function matrix obtained above and detailed in Equation (81) for the MIMO system under analysis. Also, the fact is that the system may be of a huge size since the partial fractions of the transfer functions obtained by modal series are the result of the combination of all the possible second order modal combinations. Considering these modal combinations, the poles can originate a transfer function of almost a 6750 order in the system under study. Of course, a higher number is very complicated to handle by a Bode and Nyquist routine. In order to obtain an easier expression, a simple analysis based on the highest residues was applied to the nonlinear terms. Figure 12 shows the residues location in the complex plane. Thus, due to the excessive number of residues and therefore, the same number of poles, only those with highest absolute magnitude were considered to characterize the nonlinear transfer function. The zoomed Figure 12 details the residues that were out of bounds, that is, residues whose magnitude were 0.1 i r < .
Agreeing to the reasoning of residues, the corresponding poles (modal combination) referred to the resultant residues are shown in Figure 13. This analysis could be compared with the proposal of dominant poles calculation [20,25] where the main residues and their corresponding poles are obtained using an iterative process. With respect to the complete poles of the transfer functions, the main frequencies detected in the wide range spectrum were 4. 16  (1,2) (s) for the case study of the test power system 3SM-9BUSES.
In order to obtain an easier expression, a simple analysis based on the highest residues was applied to the nonlinear terms. Figure 12 shows the residues location in the complex plane. Thus, due to the excessive number of residues and therefore, the same number of poles, only those with highest absolute magnitude were considered to characterize the nonlinear transfer function. The zoomed Figure 12 details the residues that were out of bounds, that is, residues whose magnitude were |r i | < 0.1.
Agreeing to the reasoning of residues, the corresponding poles (modal combination) referred to the resultant residues are shown in Figure 13. This analysis could be compared with the proposal of dominant poles calculation [20,25] where the main residues and their corresponding poles are obtained using an iterative process. With respect to the complete poles of the transfer functions, the main frequencies detected in the wide range spectrum were 4. 16       A spectrum characteristic of the nonlinear transfer function was obtained; the total frequency content of the nonlinear transfer function of order 6750th was analyzed, obtaining the most dominant content of some of the terms, i.e., through order 300. It was clear that there was a predominant frequency detected in the spectrum (2.08 Hz) which was a low frequency due to the electromechanical modes of the power system. Once the reduction in the nonlinear transfer function order was performed, an approximate frequency content was obtained, which is shown in Figure 14. From this figure, the similarity on frequency content with respect to the full spectrum already commented was observed. The reduced model kept almost the same frequency characteristics of the full system. Finally, Figure 15 shows the Bode and Nyquist graphs of some of the main nonlinear transfer functions, where the frequency content of each of them can be observed. In some cases, for instance ( ) ( ) Finally, Figure 15 shows the Bode and Nyquist graphs of some of the main nonlinear transfer functions, where the frequency content of each of them can be observed. In some cases, for instance Ω 2 (1,3) (s), various resonant frequencies are presented, being the predominant frequencies in the range of 2 to 3 Hz. Finally, Figure 15 shows the Bode and Nyquist graphs of some of the main nonlinear transfer functions, where the frequency content of each of them can be observed. In some cases, for instance ( ) ( )

Conclusions
In this contribution, the theory related to the linear transfer functions concepts and its extension to the transfer function with nonlinear characteristics through the application of the modal series method, as the basis of nonlinear system expansion, assuming an input force response was described in detail.
The proposed methodology has the great advantage of analyzing in both time and frequency domains, in a numerical and analytical way, the transfer functions with nonlinear characteristics of the power systems. Application of the multidimensional Laplace transform and association of variables techniques are the core of the modal series technique, which allows a closed form analytical solution of the nonlinear system, with the extension of including the control forced response. In addition, theory based on Volterra series expansion was applied in order to determine in a systematic framework the linear and nonlinear contributions of the power system model when an input signal was applied.
In the paper, the transfer function residues of transfer functions were calculated in an analytical way. Thus, the proposed methodology was successfully applied to the synchronous machine-infinite busbar classical model and to the three synchronous machine-nine buses test power system. The first case was a kind of single input-single output model; Bode analysis and Nyquist diagrams were performed in order to obtain the frequency characteristics and stability properties for the linear and nonlinear characteristics transfer functions. In addition, the time domain analysis was performed when a step input response was applied remarking the frequency oscillations and varying damping conditions presented in a more realistic operation. The second study case widely discussed along the paper was a fourth order nonlinear model, a multi input-multi output system (MIMO), in which the reference power mechanical output and voltage regulator reference were taken as the input variables. Again, the Bode analysis and Nyquist diagrams were obtained, and in addition the most dominant poles presented in the transfer functions through the pole spectrum were analyzed.
Despite the size of the nonlinear analysis which tends to be of a higher order, the contribution of the paper was oriented to the advantages of applying an analytical technique which allows several physical characteristics in both linear and nonlinear parts to be obtained. It is clear that when dealing with larger scale power systems, the complexity of analysis is increased, with necessary inclusion of sparsity techniques, effective algorithms of eigenvalues determination, etc. Future developments based on this contribution will take this problem adding the determination of dominant poles and the analysis of larger-scale power systems.

Conclusions
In this contribution, the theory related to the linear transfer functions concepts and its extension to the transfer function with nonlinear characteristics through the application of the modal series method, as the basis of nonlinear system expansion, assuming an input force response was described in detail.
The proposed methodology has the great advantage of analyzing in both time and frequency domains, in a numerical and analytical way, the transfer functions with nonlinear characteristics of the power systems. Application of the multidimensional Laplace transform and association of variables techniques are the core of the modal series technique, which allows a closed form analytical solution of the nonlinear system, with the extension of including the control forced response. In addition, theory based on Volterra series expansion was applied in order to determine in a systematic framework the linear and nonlinear contributions of the power system model when an input signal was applied.
In the paper, the transfer function residues of transfer functions were calculated in an analytical way. Thus, the proposed methodology was successfully applied to the synchronous machine-infinite busbar classical model and to the three synchronous machine-nine buses test power system. The first case was a kind of single input-single output model; Bode analysis and Nyquist diagrams were performed in order to obtain the frequency characteristics and stability properties for the linear and nonlinear characteristics transfer functions. In addition, the time domain analysis was performed when a step input response was applied remarking the frequency oscillations and varying damping conditions presented in a more realistic operation. The second study case widely discussed along the paper was a fourth order nonlinear model, a multi input-multi output system (MIMO), in which the reference power mechanical output and voltage regulator reference were taken as the input variables. Again, the Bode analysis and Nyquist diagrams were obtained, and in addition the most dominant poles presented in the transfer functions through the pole spectrum were analyzed.
Despite the size of the nonlinear analysis which tends to be of a higher order, the contribution of the paper was oriented to the advantages of applying an analytical technique which allows several physical characteristics in both linear and nonlinear parts to be obtained. It is clear that when dealing with larger scale power systems, the complexity of analysis is increased, with necessary inclusion of sparsity techniques, effective algorithms of eigenvalues determination, etc. Future developments based on this contribution will take this problem adding the determination of dominant poles and the analysis of larger-scale power systems.