Probabilistic small-signal stability analysis of power systems based on Hermite polynomial approximation

This paper proposes a probabilistic small-signal stability analysis method based on the polynomial approximation approach. Since the correct determination of unknown coefficients has a direct effect on the accuracy of the polynomial approximation method, this paper presents a method for determining these coefficients, with high coverage on the probabilistic input domain of the problem. In this method, by increasing the number of random input variables, the proposed method can continue to maintain its efficiency. After determining the unknown coefficients, the load flow results and the system state matrix are determined for random changes of all loads based on the Hermite polynomial approximation. Then, the small-signal stability of the system is probabilistically evaluated based on stochastic analysis of the system eigenvalues. The consistency and validity of the proposed method are demonstrated based on the simulation studies in the MATLAB® software environment. In the simulation studies, the performance of the proposed method is examined by comparison with Point Estimation, Cumulant, Monte Carlo, and Chebyshev polynomial-based methods, for IEEE 14-bus and IEEE 39-bus test systems. Probabilistic small-signal stability analysis of power systems Modeling of governing equations of power system based on Hermite polynomial approximation Forming the Collocation matrix based on a robust method Probabilistic small-signal stability analysis of power systems Modeling of governing equations of power system based on Hermite polynomial approximation Forming the Collocation matrix based on a robust method


Introduction
The increasing expansion of the power systems has caused it to be affected by more random factors. The more complex the system, the increase in load changes, the increasing penetration of renewable resources, the development of electric vehicles, the increase in nonlinear loads, etc.
are among the mentioned influencing factors. Stochastic changes in power load can be a factor threatening the stability of the power system. These changes are continuous and their impact on sustainability is of great importance in the control, operation, and planning of power systems. small-signal stability analysis is a fast and efficient method for evaluating the condition of power systems based on the position analysis of eigenvalues. Since the power system is a multivariable nonlinear system, to perform smallsignal stability analysis, the system state equations are linearized around the operating point and based on this, the system state matrix is formed and its eigenvalues are calculated and taken to analysis. However, because the power consumption often has uncertainty and random changes, the operating point of the system also shifts randomly due to load changes. Therefore, the use of conventional methods in which system eigenvalues are determined and analyzed based on a specific operating point, cannot lead to a comprehensive and reliable assessment.
In contrast, if the behavior of the system is predicted and analyzed probabilistically with respect to the effects of random changes of uncertain parameters, it can provide a more accurate assessment of the state of the power system and prevent technical and economic decisions from being too cautious and non-optimal [1]. So far, various methods have been used to study random phenomena such as robust optimization method, interval-based method, feasibility study method, and probabilistic method for modeling random changes, the most effective of which are probabilistic methods [2,3]. In a probabilistic method, the random changes of the input variables are modeled as the probability distribution function and based on that, the cumulative distribution function (CDF) or the probability density function (PDF) of the random output variables are calculated [3].
The probabilistic methods presented so far can be classified into two main categories of numerical and analytical methods. One of the common numerical methods for probabilistic analysis of stochastic processes is the Monte Carlo simulation (MCS) method. In this method, a large number of points of random inputs are selected and for each of them, the stability problem is definitively investigated. Then, based on the results, PDF functions, CDF functions, or other measurement criteria such as average or standard deviation of random output variables are calculated [4]. Due to the very high computational burden and time-consuming feature of the Monte Carlo method, this method is currently used as a criterion for checking the accuracy and precision of other probabilistic analysis methods. Recently, improved numerical methods such as sequential and quasi-sequential Monte Carlo methods have been introduced, which have slightly reduced the computational burden, but their use is still very time consuming [5].
One of the analytical methods to study the small-signal stability of power systems is the Point Estimation method (PEM) [6,7]. In PEM, load flow calculations are performed for a certain number of operating points, and based on that, the average, standard deviation, and probabilistic moments of the output variables are calculated. Although the PEM can maintain the nonlinearity of system equations, its accuracy in estimating high-order moments of the probability distribution function of output variables is not very high, especially for systems that are complex and have large input random variables [7].
Cumulant method (CM) is another common method of probabilistic analysis in power systems in which a linear model between output and input random variables is determined by linearizing the system equations around the nominal operating point [8,9]. Then, using the cumulants of input random variables, and the linear model of the system, the cumulants of random output variables are calculated. The most important advantage of the cumulant method is that it has a very small calculation burden and therefore has a high execution speed. However, due to the linearization of the equations governing the system around the nominal operating point, when the input random variables have significant changes to their nominal value, the accuracy of this method is affected [10].
In addition, since it is not possible to directly calculate the probability distribution functions of the output random variables in PEM and CM, to approximate the PDF and CDF functions of the variables, series such as Gram-Charlier, Cornish-Fisher, or Edgeworth series should be used [6,9]. However, studies show that the probability distribution functions obtained from these methods are not accurate and even in some cases, the resulting CDF function produces negative values [11].
In [12], since the analytical determination of the smallsignal stability region is very difficult, a hyper-plane-based approximation method is presented to depict the robust small-signal stability boundary of a power system. In [13], machine learning-based algorithms are used for the probabilistic assessment of dynamic security in a power system. In [14], to probabilistic analysis of small-signal stability in a power system, a framework is presented based on Gaussian learning.
The polynomial approximation is a method that has recently been proposed to overcome the shortcomings of numerical and analytical methods in the study of random phenomena in the power systems and has yielded acceptable results. In [15] and [16], the probabilistic load flow calculations of a distribution network under conditions of random load changes have been performed based on polynomial approximation, and the efficiency of this method has been investigated by comparison with the Monte Carlo method. The polynomial approximation method has also been used in the planning [17] and dynamic analysis [18] of power systems, which has acceptable results and good computational speed.
In the polynomial approximation method, the nonlinear equations governing the power system are approximated using a chaotic polynomial. The input variables of the mentioned polynomial are input random variables of the power system, and its output variables depending on the type of problem, are random variables, which their probabilistic properties such as average, standard deviation, or probability distribution functions are among the objectives of the problem in the power system study. The most important advantage of this method is that without using numerical series, the probability distribution functions of random output variables can be directly obtained, which in addition to high speed have good accuracy. However, the proper performance of the polynomial approximation method depends on several key factors, the most important of which is the selection of the type and order of the polynomial, as well as the proper determination of its unknown coefficients [19]. There are several polynomials that can model the governing equations of probabilistic systems. However, among them, the Hermite polynomial is more in line with the conditions of power systems, and the range of its random variables changes [19].
In [20], to investigate the effect of random power changes of a solar power plant on the power system's small-signal stability, the polynomial approximation method has been used. However, since the polynomial proposed in [20] is a first-order Hermite polynomial, the results are not very accurate. Studies in [21][22][23][24] have shown that increasing the order of the polynomial, in addition to accuracy, also greatly increases the burden of calculations. For a reasonable burden of calculations, the use of second-order Hermite polynomials shows acceptable efficiency and accuracy [21][22][23][24].
Another important issue in the polynomial approximation method that has a significant effect on its efficiency is the determination of its unknown coefficients. In [21], a second-order Hermite polynomial is used to investigate the voltage stability of a power system. In order to determine the unknown coefficients of polynomials in [21], the Galerkin method has been used. In the Galerkin method, the nonlinear equations governing the system and the polynomials that approximate the equations are combined to form a set of polynomial equations. By solving this system of equations, the desired unknown coefficients are determined. However, determining unknown coefficients using the Galerkin method has two major drawbacks. First, the use of this method causes a change in the equations of the problem and therefore it is not a comprehensive solution and for each problem, equations of the system must be formed and solved separately. The second drawback is that, as the dimensions of the problem increase and the number of random input and output variables increases, the complexity of the equations increases exponentially, increasing the volume of calculations and the time required to solve the problem.
In [22], a polynomial approximation method is presented for probabilistic small-signal analysis in power systems. However, while there is no mention of how to determine the unknown coefficients of the polynomial, the computational accuracy and efficiency of this method are only demonstrated for two uncertain input variable.
In [23], the probabilistic load flow problem of power systems is approximated using second-order Hermite polynomials and its unknown coefficients are determined using the Collocation method. In the Collocation method, using the roots of Hermite polynomial with a higher order, certain points of the input random variables are determined, and based on them, the equations of the problem are deterministically solved. In other words, using this method, the inputs and outputs of Hermite polynomials are determined for several specific points that can be used to calculate unknown coefficients.
In [24], the probabilistic analysis of the small-signal stability is performed based on the polynomial approximation method. The Hermite polynomial coefficients used in [24] are also determined based on the same method and similar to the process used in [23]. This method, unlike the Galerkin method, does not change the equations of the system and is, therefore, a comprehensive method. However, the problem with the Collocation method process in [23] and [24] is that it only works for a limited number of random input variables, and as the number of variables increases, the time and burden of calculations increase exponentially, and in some cases even it leads to divergence.
In this paper, the probabilistic small-signal stability of a power system is analyzed based on the polynomial approximation method. In the proposed method, the governing equations of power systems are approximated based on second-order Hermite polynomials, and the system eigenvalues are calculated using the approximated load flow results as well as the system state equations. The main contribution of this paper is to present a novel approach based on the Collocation method for determining the unknown coefficients of the Hermite polynomials. The proposed approach in this paper can determine polynomial coefficients without divergence for any number of random input variables. In this method, with the burden and time of calculations much less than previous methods and without using approximations by series, the probability distribution functions of system eigenvalues have been determined for random changes of power consumption in all loads and resulting operating point displacements. The consistent and effectiveness of the proposed method is evaluated in contrast with a numerical (i.e. Monte Carlo) method, two analytical (i.e. CM and PEM) methods, a Chebyshev polynomial-based approximation method, and the methods presented in [23,24].
The outline of the paper is as follows. Section 2 presents the modeling of synchronous generators, polynomials approximating of governing system equations, and | https://doi.org/10.1007/s42452-021-04765-4 probabilistic small-signal stability analyzing are presented. Based on the proposed method, the probabilistic smallsignal stability of IEEE 14-bus system, as well as IEEE 39-bus system are analyzed and discussed in Sect. 3. Finally, the concluding remarks are presented in Sect. 4.

Reduced-order model of synchronous generator
When the power system is in steady-state conditions, its state equations can be linearized around the operating point and the system eigenvalues can be determined based on the obtained equations. In this paper, the smallsignal stability of the power system is analyzed based on the reduced order 3 model of the synchronous generator, together with the excitation system model, for two reasons: 1) the critical eigenvalues obtained from the 3rdorder model of the synchronous generator are not much different from the values obtained from the 8th-order model [25], 2) what is approximated by the polynomial method is the probabilistic load flow calculations necessary to determine the operating point of the system, as well as the state matrix, therefore, the order of synchronous generator model has no noticeable effects on the accuracy of the proposed method. The mentioned model is expressed by equations (1)-(4) [25].
where, is the rotor angle, is the rotor speed, E ′ q is the transient voltage along the quadrature axis, E fd is the excitation voltage, V t is the magnitude of the terminal voltage, V ref is the reference value of the terminal voltage, and I q and I d are the generator current along with the quadrature and direct axes. T ′ d0 is also the time constant of the direct axis of the synchronous generator in the transient state, and T A and K A are also the time constants and the gain of the excitation system, respectively. By rewriting equations (1)-(4) and organizing them into a matrix form, we will have [25]: In these equations, the coefficients K 1 to K 6 and also Δ are: Now, the eigenvalues of each of the generators can be calculated using the state matrix A and the same matrix I and based on the relation det( − IA) = 0 , where, = ± j are system eigenvalues.

Approximation of system equations based on Hermite polynomials
In this paper, the power consumption of all loads of the system is considered uncertain and, as usual, is modeled with the normal probability distribution function [26]. In this modeling, the mean value of the i th load ( i ) is considered equal to its nominal value and the standard deviation of the i th load ( i ) is considered equal to 10% of its mean value [10]. To use the polynomial approximation method, the non-deterministic inputs of the problem need to be converted to standard normal random variables [19]. Therefore, in order to convert the load changes from the normal distribution to the standard normal distribution and vice versa, equations (16) and (17) have been used: where, x i is the i th random input variable of the problem, ( i th load) and i is a standard normal random variable. If Y is the vector of the load flow results, the load flow equations can be approximated using the Hermite polynomial, as follows [21,24]: where, m is the number of random input variables and k j is the unknown polynomial coefficients. These coefficients must be determined in such a way that the above equation can determine the load flow results for the values of the input random variables (load changes). In (18), Γ n is a Hermite polynomial of order n, the general equation of which is determined by (19) [21,24].
For example, Hermite polynomials from orders n = 0 to n = 3 can be expressed by the following equations: Comparing the accuracy of different orders of Hermite polynomials has shown that the response obtained from second-order polynomials is slightly different from the response obtained from Hermite polynomials with the order of 3 and above [21,24]. Therefore, in this paper, second-order Hermite polynomials have been used to model the load flow calculations. If Y l is the l th output variable of the load flow problem, it can be expressed in terms of input variables as:

Determining the unknown coefficients of Hermite polynomials
To calculate the unknown coefficients k 0 , k j , k jj , and k ij in (21), it is first necessary to determine some sample points, with certain values of . In an n-order Hermite polynomial, sample values should be randomly selected from a permutation of zero and roots of (n+1)-order Hermite polynomial. The number of sample points also depends on the order of the polynomial and the number of random variables of the system and is determined by (22) [27].
where, r, n, and m are the number of sample points, the order of Hermite polynomials, and the number of random variables, respectively. By rewriting Equation (21), a matrix equation can be formed as follows: where, [K] is the vector of unknown coefficients and [Y] is the vector of load flow results. The l th row of (23) can be expressed for m input random variables as follows: As can be seen in (24), the matrix [M] is an (r × r) square matrix whose l th line is the combination of the number of one and values , which is one of the operating points of the system, is called the Collocation point [27]. It should be noted that to determine the unknown coefficients [K], it is necessary to select the r Collocation point in such a way that after forming the matrix [M] and replacing it in (23), create a set of independent linear equations, so that the problem has only one unique answer. In other words, the Collocation points must be selected so that the rank of the matrix [M] is complete and equal to r. Since the accuracy of load flow calculations using polynomial approximation is highly dependent on the correct determination of unknown coefficients, it is necessary to select the Collocation points in such a way that in addition to the completeness of the matrix rank of [M] , have maximum coverage on the m-dimensional probabilistic input space of the problem. The method presented in [23,24] is that first a large group of Collocation points is formed and then by  [23,24] has been studied only for three random input variables and has acceptable results. However, this method only operates when the number of input variables in the problem is small and r is not a large number. Because when r increases, the dimensions of the matrix [M] increase by the same amount. Since its rank must be calculated in each iteration of the matrix formation process, it will be very time-consuming to use the method presented in [23,24]. In addition, if the iterative process is complete and a full rank matrix is not yet formed, the method presented in [23,24] will practically diverge.
In this paper, a new method is proposed for forming the matrix [M] and calculating the unknown coefficients of Hermite polynomials. In the proposed method, first, a group of Collocation points is formed with a random combination of zero numbers and roots of third-order Hermite polynomial. A row of [M] matrices is then formed for each member of the group. Since it contains R ≫ r , by juxtaposing these rows, the matrix [M g ] with dimensions of (R × r) is obtained. Then the rank of the matrix [M g ] is calculated and if its rank is r, by deleting its non-independent rows, the matrix [M] with the full rank is formed and the process ends. To ensure convergence, if the rank of the matrix [M g ] is less than r, a new group of Collocation points is formed and this process is repeated. Although determining the rank of the matrix [M g ] has more calculations than determining the rank of the matrix [M], it should be noted that the rank determination calculations in the proposed method only once done. The method proposed in this paper, while simple and feasible, has a high computational speed, its efficiency does not decrease due to increasing the number of random input variables, and in addition, it is practically impossible to diverge.
The elements of the vector [Y] in the system of (23) are obtained by performing load flow calculations based on deterministic methods. In the above calculations, the Collocation points that form the matrix [M] are first converted into real variables (load power) using (16), and then, based on that, the load flow results are obtained. The flow chat of this procedure is shown in Fig. 1. It should be noted that only for the r number of the sample operating points, the deterministic load flow calculation is done. Then, the unknown coefficients are calculated based on the following equation.

Probabilistic small-signal stability analysis
Once the unknown coefficients [K] are determined in (23), the equation can be used to repeatedly calculate the probabilistic small-signal stability results for random changes of loads instead of the usual long calculations. For this purpose, first, the random values of the load are converted to standard normal random variables using (17) and then the load flow results are determined based on (24). Now if the system state matrix is calculated using the load flow results, the system eigenvalues can be determined based on det( I − A) = 0 . In the following, it will be shown that the use of the proposed method causes the probabilistic analysis of the small-signal stability of power systems to be performed with appropriate accuracy and very high speed.

Results and discussion
In this section, the probabilistic small-signal stability of IEEE 14-bus system, as well as IEEE 39-bus system are analyzed using the method proposed in this paper. The analysis is based on coding in the MATLAB software environment and performed using a PC with Intel ® Core ™ i7 CPU at 2.5 GHz, 8 GB RAM, and 64-bit Windows 10 OS. The single line diagram of the power systems under study is shown in Fig. 2. In order to verify the accuracy of the proposed method, the probabilistic calculations of the studied systems have been performed using the Hermite polynomial approximation method, and its results (denoted as "HAM") is compared with the analytical methods of Point Estimation (denoted as "PEM") and Cumulant (denoted as "CM"), and the numerical method of Monte Carlo simulation (denoted as "MCS"). In order to further increase the accuracy of the MC method, as a criterion method, 10,000 and 30,000 points have been selected from the probability distribution function of buses load powers in the 14-bus and 39-bus systems, respectively. In addition, the computational efficiency and convergence of the proposed method have been evaluated in comparison with the method presented in [23,24], as well as the Chebyshev polynomial approximation method (denoted as "ChAM").

IEEE 14-bus test system
In Table 1, the mean and standard deviation values of the voltage magnitude of busbars obtained from the proposed method are compared with the MCS, PEM, CM, and ChAM methods. Since bus 1 is of slack type and buses 2, 3, 6, and 8 are also of PV (voltage-controlled) type, their voltage magnitude did not change in the face of random changes of load powers, and therefore, their standard deviation in all methods is zero. As can be seen, the average values of the voltage magnitudes obtained from the polynomial approximation methods (HAM and ChAM), for all PQ buses of the system, have a very small difference (in the range of 10 −5 pu) with the values obtained from the MCS method. In addition, comparing the standard deviation of the HAM and ChAM methods with MCS method confirms that the polynomial approximation methods have acceptable accuracy. It can also be seen from Table 1 that the voltage magnitude of bus 14, |V 14 | , has the highest standard deviation, and the voltage magnitude of bus 12, |V 12 | , has the lowest standard deviation among all busbars. This indicates that the load uncertainties have had the greatest impact on |V 14 | and the least impact on |V 12 | . The cumulative probability function of the voltage magnitude of busbars 14 and 12 for random changes of all loads of the system is shown in Figs. 3(a) and (b), respectively. As can be seen, the results of the approximation methods of HAM and ChAM, unlike the analytical methods of PEM and CM, are slightly different from the results of the Monte Carlo method. It is noteworthy that the average voltage of the busbars in all five methods is almost the same, and what has caused the difference between the results is the significant difference between the standard deviation of the PEM and CM methods with the MSC and HAM methods. Although the results of the ChAM method are more accurate than the analytical methods, they have lower accuracy than the MCS and HAM methods. Therefore, it can be said that the method proposed in this paper, in addition to being faster, has also shown good accuracy in calculating the probabilistic load flow of the power system.  Table 2 shows the mean and standard deviation of the eigenvalues of the system under study. As can be seen, the mean values obtained from the proposed method are not significantly different from the results of other methods and are almost equal to each other. Examination of standard deviation values shows that the results of the proposed method have an acceptable agreement with the results of the MCS method. However, the standard deviations obtained from the PEM and CM are slightly different from the other MCS and HAM methods. The reason for this phenomenon is the estimation of the cumulative probability function in the PEM and CM based on the Cornish-Fisher and Gram-Charlier series, respectively. Although the PEM and CM are the best among analytical methods, the proposed method offers a much better estimate. Figure 4 shows a graph of the mean eigenvalues of the system under study on the plane of the complex numbers. As expected, according to the dimensions of the state matrix, the system eigenvalues can be clustered in four modes. The first mode is completely real eigenvalues that are close to zero. The second and third modes are conjugated to each other and are considered critical eigenvalues of the system due to their damping factors. The eigenvalues of the fourth mode are also completely real, which are more dependent on the excitation system and are the farthest eigenvalues from the imaginary axis. In modes 2 and 3, the eigenvalues for generators 1, 2, and 3 have the lowest damping factors. Re-examination of Table 2 shows that among these three generators, generator 3 has the lowest damping factors. It is noteworthy that Generator 3 also showed the highest standard deviation. Therefore, in addition to the lowest stability margin, generator 3 also has the highest sensitivity to load power changes among system generators. Figure 5(a) shows the cumulative distribution function of the real part of the mode#2 eigenvalue of the third generator ( 2,G3 ) . As can be seen, the results of the proposed method are not significantly different from the MCS method and have good accuracy. However, the CDF function obtained from the PEM and CM methods, although have good accuracy in terms of the mean value, has an asymmetric difference in other points with the results of the proposed and MCS methods. In fact, it should be noted that based on the results of the HAM, ChAM, and MCS methods, the CDF function of the real part 2,G3 has asymmetry or skewness from the right, however, the PEM and CM methods have not been able to approximate this skewness correctly. Figure 5(b) shows the probability distribution function of the real part 2,G3 . As can be seen, the PDF function shown in Fig. 5(b) has positive skewness in the HAM, ChAM, and MCS methods. The value of this skewness is 0.4421 but the skewness value obtained from the PEM and CM methods are 0.4079 and 0.4082, respectively, which are significantly different from the correct amount. The functions of cumulative distribution and probability distribution of the real part of the mode#1 eigenvalue of the third generator 1,G3 are also shown in Figs. 6(a) and (b), respectively. As can be seen, the results of the HAM, ChAM, and MCS methods are not significantly different, but the PEM and CM methods still failed to correctly approximate the CDF and PDF functions of the real part 1,G3 . As can be seen in Figs. 6(a) and (b), the CDF and PDF functions approximated by the HAM, ChAM, and MCS methods have a negative skew (left), however, in the PEM and CM methods, this skew has not been approximated correctly. Another noteworthy point in Figs. 3, 5, and 6 are that the probability functions obtained from the proposed method are also completely consistent with the MCS method in terms of kurtosis, but the results obtained from the PEM and CM methods are almost erroneous from this point of view. By comparing the ChAM and HAM methods, it can be seen that the accuracy of the proposed method is higher and the results are more consistent with the MCS method.

IEEE 39-bus test system
The graph of the mean eigenvalues of the IEEE 39-bus system on the plane of the complex numbers is shown in Fig. 7. In this system, similar to the 14-bus system, the system eigenvalues can be clustered in four modes. The first and the fourth modes are completely real eigenvalues that are close to zero and far from the imaginary axis, respectively. The second and third modes are conjugated to each other and are considered critical eigenvalues of the system due to their damping factors. In modes#2 and #3, the eigenvalues for generators 9, 6, 2, 7, and 3 have the lowest damping factors.
The mean and standard deviation of the eigenvalues of the IEEE 39-bus system in modes#2 and #3 are shown in Table 3. As can be seen, the mean values obtained from the proposed method are not significantly different from the results of other methods and are almost equal to each other. Examination of standard deviation values shows that the results of the proposed method have an acceptable agreement with the results of the MCS method. However, the standard deviations obtained from the PEM and CM are slightly different from the other MCS and HAM methods. The reason for this phenomenon is the estimation of the cumulative probability function in the PEM and CM based on the Cornish− Fisher and Gram-Charlier series, respectively. Although the PEM and CM are the best among analytical methods, the proposed method offers a much better estimate.
Re-examination of Table 3 shows that among these generators, generators of 9 and 6 have the lowest damping factors. It is noteworthy that Generators 9 and 7 also showed the highest standard deviation. Therefore, in addition to the lowest stability margin, generator 9 also has the highest sensitivity to load power changes among system generators.
Figs. 8(a) and 8(b) show the CDF and the PDF functions of the real part of the mode#2 eigenvalue of the 6th generator ( 2,G6 ) . As can be seen, the proposed method exhibits an acceptable accuracy and is in consistent with the MCS method. However, the CDF and PDF functions obtained from the ChAM, PEM, and CM methods, despite good accuracy in terms of the mean value, have an asymmetric difference in other points with the results of the proposed and MCS methods. It should be noted that in 39-bus system, a fundamental difference has been shown in the accuracy of the approximation of HAM method with ChAM method.
The CDF and PDF functions of the real part of the mode#2 eigenvalue of the 9th generator 2,G9 are shown   , respectively. As can be seen, the results of the HAM and MCS methods are not significantly different, and the shape of the CDF and two peaks PDF obtained from MCS method has almost been fitted by the proposed method. However, the PEM and CM methods failed to correctly approximate the CDF and PDF functions of the real part 2,G9 . In addition, Figure 9 illustrates that as the scale and complexity of the system increase, the ChAM method is unable to approximate the governing equations correctly and loses much of its accuracy. Investigating the results of the proposed method and comparing them with the Monte Carlo method as a criterion, shows that the proposed method has considerable accuracy in analyzing the small-signal stability of power systems. The important point is that the proposed method has not only been able to accurately estimate the mean and standard  deviation of the output variables, but also the CDF and PDF functions of probability variables have approximated with high accuracy, so that even the qualitative characteristics of these functions, such as skewness and kurtosis, are also properly approximated.

Evaluation of computational efficiency
In Table 4, the computational time as well as the number of deterministic load flow calculations of the proposed method is compared with the PEM, CM, MCS, ChAM methods, and also the method presented in [23,24]. As can be seen, the number of deterministic calculations in the Monte Carlo method should be high in order for the response accuracy to be acceptable, which leads to heavy and time-consuming calculations. However, in analytical and polynomial approximation methods, since the deterministic calculations are very limited and are performed only for a certain number of operating points, the computational burden and computation time is much less than the Monte Carlo method. By comparing analytical methods (PEM and CM) and polynomial approximation methods (HAM and ChAM), it can be said that since the   number of deterministic calculations required to determine polynomial coefficients in the approximation methods is slightly more than analytical methods, the time and volume of calculations in the approximation method is slightly more than the analytical methods. However, it should be noted that the accuracy of the calculations in the approximation methods is very high and almost equal to the Monte Carlo method, but the burden of calculations and its complexity is about the analytical methods, and this feature is the biggest advantage of polynomial approximation methods over other methods. Comparing methods HAM and ChAM, it can be seen that the number of deterministic calculations and computation time in method HAM is more than method ChAM, however, the simulation results have shown that the accuracy of method HAM is also more than method ChAM. It is worth noting that the process of forming matrix [M] for the random changes of all system loads using the methods presented in [23,24] is unstable and does not result even after over ten minutes.

Conclusion
This paper proposes a polynomial approximation-based method for probabilistic small-signal stability analysis of power systems. In the presented method, the governing equations of power systems are modeled based on Hermite polynomials so that instead of the governing equations, the obtained polynomials are used to calculate the state matrix and eigenvalues of the system. The unknown coefficients of polynomials, which have a direct effect on the accuracy of this method, are determined based on a novel approach that in addition to high efficiency, also has a good speed and convergence. The proposed approach makes it possible to use the polynomial approximation method for power systems with any number of random input variables. In order to evaluate the accuracy and efficiency of the proposed method, the IEEE 14-bus and IEEE 39-bus systems have been simulated in MATLAB software and all its loads have been modeled as random variables. The validity and effectiveness of the proposed method is demonstrated in contrast with the Monte Carlo, Point Estimation, Cumulant, and Chebyshev polynomial approximation methods. Evaluating the simulation results shows that the proposed method has not only been able to accurately estimate the mean value, the standard deviation, and the probability distribution functions of the output variables, but also the properties of the mentioned functions, including kurtosis and skewness, have properly approximated. Ease of use, fewer calculations, high analysis speed, accuracy, and efficiency are some of the advantages that the proposed method has shown in probabilistic small-signal stability analysis of power systems, especially when the number of input random variables is significant.

Declarations
Conflict of interest The authors declared that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.