CC BY
8
FEASIBILITY OF OPERATING CONDITIONS AND THE LOCATION OF SENSOR VARIABLES IN THE ELECTRIC POWER SYSTEM
Irina I. Golub, Oleg N. Voitov, Evgeny V. Boloev1 •
Energy Systems Institute of SB RAS, Irkutsk, Russian Federation
e-mail: [email protected]
ABSTRACT
In the paper the methods of probabilistic load flow, including linear method of generalized disturbance, are used to detect sensor variables in an electric power system, identify their probabilistic characteristics, detect critical variables, for which the probability that they lie in the feasible range is lower than the required one, and select the control actions to increase this probability. The controls are selected by the method similar to the method of deterministic equivalent by subsequently and iteratively solving deterministic and probabilistic problems. The method of contribution factors makes it possible to choose from a set of possible controls a vector of control that contains the minimum number of components. The presented numerical results on the example of test and real networks demonstrate the efficiency of the proposed approaches to the electric power system operation control.
1. INTRODUCTION
State variables of the electric power system should lie within certain feasible limits, which meet reliability and quality requirements. The probability that the variable will go beyond its feasible limits in case of disturbances depends on the sensitivity of the variable to disturbances, a feasible range of its variation, and proximity of its current value to the limit.
Such network elements whose state variables change largely due to random external disturbances are called sensors [1]. Inhomogeneity of the electric power system that leads to the emergence of sensors is determined by both the operating conditions of the system and parameters of the network elements that are called weak places in [2].
Two approaches can be used to detect sensors and weak places. The first approach is related to the analysis of disturbance scenarios and responses of variables to the disturbances that are expressed, for example, by deviations of variables.
In the first approach the sensitivity of a variable to external disturbances can be determined on the basis of singular or spectral decomposition of the Jacobian matrix. If the spread of singular values (eigenvalues) is large, it indicates severe inhomogeneity of the electric network. The greater the difference between the first singular value of the Jacobian matrix and the rest of the values, the more reasons to make a conclusion about the behavior of the variables on the basis of the first singular decomposition summand connected with the minimum singular value.
1 The study has been supported by grant 1507.2012.8 Leading Scientific School.
The idea of using the maximum components of the eigenvector that correspond to the minimum eigenvalue to choose the most informative set of measurements [3] was applied to detection of sensor nodes in an electric network by means of the spectral analysis in [1], [4] and singular analysis in [5].
The variables most sensitive to external disturbances can be identified by using a scalar value of the first generalized disturbance. The greatest contribution to the first generalized disturbance is made by the disturbances at nodes corresponding to the components of the left singular vector, while the components of the first right singular vector distribute the generalized disturbance among the nodes of the electric network.
In the second approach used in this study the response of a random variable to disturbances can be determined by the methods of probabilistic load flow, including the linear analytical method which uses the scalar variance value of the generalized disturbance. The numerical characteristics of the variables obtained in the calculation allow us to find the probability that the variables are feasible.
A large response of variables to a disturbance is significant in the case if it changes some criterion of power system operation, for example the criterion of operation feasibility which is considered in this study.
To increase the probability that the variable lies in the feasible region it is necessary either to reinforce the network, which will improve the Jacobian matrix conditioning, reduce the response of the variables to the disturbance, and expand the feasible region, or to find the appropriate controls to make the variable mean value shift inside the feasible region.
The required probability is provided by iteratively and consecutively solving the problem of probabilistic load flow and the problem of determining feasible operating conditions by the deterministic equivalent method, which implies searching for a feasible solution with a shift of the critical variable mean value inside the feasible region. The critical variables are the variables for which the probability to fall within the specified interval is less than the required one. If the solution exists, a control vector with the minimum number of components is chosen from a set of possible controls, using the method of contribution factor [6], and the required increment in the control vector is determined, which increases the probability that the critical variable lies in the feasible region.
Analysis of contemporary methods for the calculation of probabilistic load flow and their use to solve various power engineering problems is presented in the overviews [7] and [8].
The methods of probabilistic load flow can be divided into the methods of linear and nonlinear approximation, which can be both iterative and noniterative, and numerical methods. The noniterative linear methods include the method of moments [9] that are formed on the basis of the Jacobian matrix, and the method of convolution [10]. The foundations of the linear iterative method that was called the method of statistical linearization were developed in [11] and successfully applied in [12]. In the noniterative [13] and iterative [8] nonlinear methods the moments are formed on the basis of the Jacobian and Hessian matrices.
The method of Monte Carlo and point methods [14], in which ordinary programs for deterministic load flow are used for the calculations of probabilistic load flow are referred to the numerical methods.
The study to be mentioned among the first to consider constraints in the calculation of probabilistic load flows is [15]. The need to solve the indicated problem by the procedure for calculating constrained optimal power flow has resulted in the development of the methods that combine deterministic and probabilistic approaches. The theoretical foundations of this approach called the method of deterministic equivalent are presented in [16].
2. DETECTION OF SENSOR VARIABLES BY THE METHODS OF PROBALISTIC LOAD FLOW
The mean square deviations of nodal voltage magnitudes and phases can be determined by the specified mean square deviations of nodal powers in the linear analytical method, using the expression relating the changes in phases A8 and magnitudes AU of nodal voltages and the changes in active A P and reactive A Q powers in the system of linear equations
, fAP
=J
VAU ,
AQ,
(1)
where J - Jacobian matrix.
The means mA5 AU and covariances m2A5AU of changes in the magnitudes and phases of
voltages are determined through the means mAP, AQ and variances of loads m2AP,AQ at the point of
solution to the nonlinear system of steady state equations of electric power system as
mA5,AU = J mAP,AQ, (2)
M2A5,AU = J M2AP,AQ(J ) . (3)
The expressions of the means and covariances of changes in the active and reactive power flows as well as differences of voltage magnitudes and phases can be written in an analogous way. In particular for covariances they will have the form
M2Ap,,AQt] = JijJ M2AP,AQ (j i jJ ) , (4)
M 2 A (5,-5 ] ),A(U,-U] ) =MTJ M 2 A P, A Q (j ^ ) ' M , (5)
where Jn - matrix of partial derivatives of active and reactive power flows in the tie ij with respect
to magnitudes and phases of nodal voltages, M - the incidence matrix.
The numerical characteristics of loads under the assumption about their normal distribution can be obtained by using the Laplace function [9], also called the error function.
The linear method of generalized disturbance is based on the combination of linear analytical method (2), (3) with the method of singular analysis which implies singular decomposition of the asymmetrical Jacobian matrix
J=WZVr =Jwj a
■ ■ vr,
(6)
where W = (w1,
j=1
Wl,W2,...,Wn
) and V =(v1,v2,...,vn ) - orthogonal matrices, whose columns represent left
and right singular vectors, and E - diagonal matrix of singular values <a2 <a3 <.
arranged in the ascending order.
Taking into consideration decomposition (6), expression (1) can be written in the form
fA P'
. <a,.
fAô> =J 1
vAU.
..T ^AP
'i v
A Q.
= & A S
<0
(7)
^ w;
where À S^^ - the i-th generalized disturbance.
If the first singular value =amin is considerably lower than the rest of the singular values
the largest contribution to the changes in phases and magnitudes of nodal voltages is made by the first term of the sum (7)
f Aô
AU
= v,
w,
a,
fA P
A Q.
= v A S
(1)
(8)
The first generalized disturbance A Sis the same for all nodes and is distributed among state variables in proportion to the components of the first right singular vector Vj. In this case the
maximum contribution to the first generalized disturbance is determined by the disturbances at the nodes that correspond to the maximum components of the first left singular vector w1.
The number of disturbance variants is infinite. They may differ both in composition and in value. The linear method of generalized disturbance does not require the scenario of change in the nodal powers to be specified and makes it possible to assess a set of disturbance scenarios by using the specified value of covariance ^ (I) of the generalized disturbance. The most significant
disturbances are the generalized disturbances that correspond to the first minimum
(1) _V1W!
M2 AS, AU M* 2
T f T\
2AP,AQ
u1
or k minimum singular values that are close to one another
_V1 M2A^VT (9)
k
^AS.AU = ZV< ^2A. (10)
i = 1
Similar expressions including the variance of generalized disturbance can be written for the numerical characteristics of changes in power flows, differences of nodal voltage magnitudes and phases.
3. PROBABILISTIC CONSTRAINED LOAD FLOW
A strong response of the sensor variable to the disturbance is not dangerous in itself if the variable remains within feasible limits after the disturbance. Therefore, the most important indicator is the required value of probability that the random value lies in the specified feasible range.
In the case that the calculation of probabilistic load flow results in critical variables for which the probability to fall within the specified interval is less than the required one, the probability can be increased either by decreasing the mean square deviation of the variable or by shifting its mean value inside the feasible region.
The mean square deviation of the critical variable being also a sensor variable, can be decreased for example by the reinforcement of weak ties [1], [4]. Another possibility is to choose the control actions that decrease the distance between the mean value and median of the distribution density curve f (x) on the feasible interval.
Such a criterion is used in the case where variable * has the law of distribution other than the normal law, and the approximated probability density curve can be obtained by three or a greater number of moments, using the Gram- Charlie expansion [9]
f(x) = £ ciH1 (x)|(x), (11)
j=0
where | (x) - probability density for normal distribution, Hj (x) - orthogonal Hermite polynomials, cj - coefficients built on the basis of the second- and higher-order moments.
The probability density curve when approximated on the basis of two moments is symmetrical. This makes it possible to transform the criterion for the selection of controls into the minimization of distance between the mean value of the variable and the center of a feasible interval of its change towards the point of the required probability value. When the constraints on the variable are specified symmetrically with respect to its nominal value the center of the interval is the nominal value of the variable.
In order to choose the controls to provide the required probability for the critical variables to lie within the feasible limits, the method similar to the method of deterministic equivalent [15], [16] is used. In this method successively deterministic and probabilistic problems are solved. Solving the deterministic problem suggests shifting the mean value to the center of the feasible interval but not
narrowing this interval for each critical variable as it is done in the method of deterministic equivalent.
The variable mean value to be obtained as a result of control can be determined by using the inverse error function. This function allows one to determine the interval As of change in the normally distributed random variable, which makes it possible at a specified value of mean square deviation to provide the required probability that this variable falls within this interval.
The interval As is compared to the known feasible interval Asfeasible of the variable change.
If 9 = (As feasible -As)> 0, then the mean of the variable should equal 9 , otherwise a conclusion about the impossibility of providing the required probability that the variable lies within the feasible interval Asfeasible and the need to shift the mean value of the variable to the center of the feasible
interval is made. In the first situation the required shift of the mean z of variable A z will equal
A z = 9 - z, and in the second situation A z = p,A z.
An algorithm for increasing the probability that the variables fall within the feasible limits is iterative, its each k -th iteration contains the following main steps.
1. The deterministic problem is solved to obtain feasible operating conditions of electric power systems subject to
W (.X,Y ) = 0, (12)
X^ < X < X_, (13)
Fmm < F(X,Y)< Fmax, (14)
Y <Y <Y (15)
mm max 5 \ s
where (12) - system of equations of nodal power balances, (13)-(15) - constraints on the dependent variables X, that include magnitudes and phases of nodal voltages and functional variables F, that contain active and reactive power flows; (15) - constraints on controls or independent variables Y, such as active and reactive power of generation and transformation ratios of tap-changing transformers.
Problem (12) to (15) is solved by combining the reduced gradient and quadratic programming methods [15], and if there exists a feasible solution for vectors Xk and Yk then the algorithm goes
to step 2. Otherwise, the algorithm stops.
2. The probabilistic load flow is calculated, the numerical characteristics of variables and probability of meeting the constraints (13)-(15) are determined. If for all the variables the required value of probability is provided, the algorithm stops. Otherwise, the number Nv of critical variables
zkj is determined for which the required probability value is not provided and an estimate of the
shift Aj of its mean which leads to an increase in the probability is calculated. If in the adjacent
A k+1-A\
< is met, where % - a set small
iterations for each critical variable zj the condition
number, the required probability values cannot be reached and the algorithm stops. Otherwise, the algorithm goes to step 3.
3. The deterministic optimization problem is solved to determine the vector of control Y„ = Yk + AYk , that provides the minimum of the criterion
mm±(z,(Y)-(zj(Yk)+A^))2, (16)
j=i
when the constraints (12)-(15) are met. If a solution to the problem is found and criterion (16) equals zero, the algorithm goes to step 4, otherwise, Yk+1 = Y„ , k = k +1, and it goes to step 2.
4. The minimum number of controls AYk is chosen on the basis of the contribution factors method [6]. This method allows the identification of the ways of transmitting active and/or reactive power from generator nodes to load nodes and the contribution of generator power to the power of flows and loads. The information about tracing the flows is used to determine the so called significant controls that affect the critical variable to the greatest extent. The significant controls underlie the formation of variants with different number of controls AYk, for each of which the solution to problem (16), (12)-(15) is searched for. When comparing the variants, the variants with the minimum number of controls are chosen. If there are several variants with equal number of controls, then the variant with the minimum active power losses is taken. Then the vector Yk+1 = Yk + AYk, k = k +1 is determined, and the algorithm goes to step 2.
3. CASE STUDY
The electric power system presented in Figure 1, which consists of 14 nodes and 15 ties, is used as a test scheme. The performance of the considered methods for this scheme is illustrated by an example of the nodal voltage magnitude control.
Figure 1. Scheme of a 14-node test network
The initial data on mean values and variances for the loads specified at all the nodes were obtained with the use of the Laplace function. Mean square deviations of nodal powers were assumed to be equal to 12 % of their mean values, which corresponds to 20 % of the load forecast error for a 0.9 probability of random value deviation from the mean.
Figure 2 presents the graphs of the mean square deviations of nodal voltage magnitudes obtained by the linear method, generalized disturbance method, and the Monte Carlo method. The graphs show that node 8 is the sensor node.
The conclusion that node 8 is the node with a sensor voltage magnitude can also be made on the basis of the singular analysis technology [5]. For this purpose the nodes can be projected on the plane in coordinates of the first and second singular vectors. The sensor nodes in such a graph will have maximum distance from the origin of coordinates. After interconnecting the nodes by the ties, the network graph projection in coordinates of the first and second right singular vectors is obtained, Figure 3.
qAU (kV ) 0102D3
1
¡8» Bsran M mn
2 4 5 6 8 100 200 202
node number
Figure 2. Mean square deviation of voltage magnitudes at the nodes of the test scheme obtained by the linear method - 1, generalized disturbance method - 2, the Monte Carlo method - 3
Figure 3. Projection of the network graph in coordinates of the first (vi) and second (V2) right singular vectors that correspond to voltage magnitudes
However, such a technology, unlike the probabilistic load flow, does not allow simultaneous identification of sensor variables and assessment of their possible variation ranges and probabilities that the variables lie within the feasible limits.
Table 1 contains the mean values and mean square deviations of voltage magnitudes obtained by the linear method of probabilistic load flow, differences between the mean values and nominal voltages, and the probabilities that voltage magnitudes fall within the feasible intervals. The feasible intervals for 500 kV voltages are taken equal to ±30 kV, and for 220 kV - ±25 kV.
Table 1. Probabilistic characteristics of voltage magnitudes at the test network nodes for the
initial state
Nodes mU <u (kV) mu - Unom (kV) P
2 522.34 3.65 22.34 0.98
4 231.49 1.32 11.49 1.00
5 512.05 6.88 12.05 0.99
6 225.17 1.86 5.17 1.00
8 508.44 17.12 8.44 0.88
100 229.24 2.03 9.24 1.00
200 528.15 4.83 28.15 0.64
202 233.62 2.05 13.62 1.00
If the required value of probability that voltage magnitudes lie within the given intervals should be not less than 0.95, the voltage magnitudes at nodes 200 and 8 can be defined as critical. Another estimate of critical voltage magnitudes A Ui can be represented by the maximum ratio of
mean square deviation qAU to the feasible variable variation range determined by the proximity of
the variable mean to its upper A Ui or lower A U. limiting value
^AU = Sa ujmin (A Ui "Mau^Mau, - AUi).
Such a possibility is illustrated in Figure 4 which shows the values of components of vector ^AU and the values of probabilities that voltage magnitudes lie within the feasible limits. Under the
initial operating conditions critical node 200 corresponds to the maximum value ^AU and minimum probability. At the nodes with the 0.98-1.0 probability that voltage magnitudes are within the feasible limits, the value of criterion ^AU does not exceed 0.5.
To provide the required probability that voltage magnitudes at nodes 8 and 200 fall within the feasible limits the two methods including the reinforcement of weak ties and selection of control actions to move the mean values of variables to the center of the feasible interval were compared. Weak ties are the ties, in which the reduction in resistances increases the minimum singular value of the Jacobian matrix, i.e. improves its conditionality and decreases the response of sensor variables to disturbances [5].
32.5 -21,510,50-
Figure 4. Values of index ^AU - 1, 2 and probability p that voltage magnitudes are within the
feasible limits - 3, 4, for the initial operating conditions 1, 3 and the conditions obtained as a result of the network reinforcement and selection of control actions 2, 4
14 12 10 8 6 4 2 0
Figure 5. Mean square deviations of voltage magnitude differences in the ties of the scheme obtained by the linear method -1,the generalized disturbance method -2, the Monte Carlo method -3
For the criterion used to detect weak ties the maximum values of the mean square deviations of changes in the voltage magnitude differences Figure 5 obtained by the linear method, the generalized disturbance method, and the Monte Carlo method are applied. According to these values, ties 5-8 and 8-200 in the test network Figure 1 are weak. The projection of the network graph Figure 3 shows that the weak ties are the longest.
+ A Ui'P H102 3 A 4 i
¡¿3 1
^s^ ^^ iii mw, %
2 4 5 6 8 100 200 202
node number
SaU, -u, ) (kV) H10 2D3J
i i
-j
m \i S ~IT~iT
\ ESTl sn S m\ s
7 T
(N CO
tie
At the initial point the mean square deviation of the voltage magnitude at node 8 amounts to Us = 17.12 kV. This means that with the probability of 0.95 the voltage will be in the interval
As = ± 33.54 kV. In the feasible interval of voltage changes equal to Asfeasible = ± 30 kV, the mean
of the variable will be 0 = (30-3354) =-3.54 < 0, whence it follows that the mean should be moved to the center of the feasible interval. This will make it possible to provide the 0.9203 probability that the voltage falls in the feasible interval.
To provide the 0.99 probability that the voltage of node 200 with mean square deviation equal to ^AU = 4.83 kV lies in the feasible interval, the centered value of the mean should equal
0 = (30-11.77) = 18.23 > 0 and the shift of the mean will be
AAU =0-^AU = 18.23-28.15=-9.92 kV.
AU200 r AU200
Since the critical variables are represented by voltages, their values are changed by the reactive power sources and transformers.
The test scheme has 11 controls with 6 reactive power sources at nodes 1, 3, 7, 101, 201, 203 and 5 tap-changing transformers in ties 2-4, 5-6, 200-201, 200-202, and 202-203.
The contribution factors method [6], [17] was used to determine the controls that make the required shifts of voltage values at nodes 8 and 200 in the test scheme by investigating the reactive power flows coming to the specified nodes, Figure 6.
Figure 6. Directions of reactive power flows arriving at nodes 8 and 200 with critical voltage
magnitudes
The analysis of the flow directions makes it possible to find the significant controls that include four sources of reactive power at nodes 1, 3, 201 and 203 and four tap-changing transformers 2-4, 200-201, 200-202 and 202-203.
The determined controls are the basis for the calculation of vector AYk with the minimum number of controls that include the reactive power of node 203 and the transformation ratios of transformers 200-201, 200-202 and 2-4.
Figure 7 presents the distribution density curves of changes in voltage magnitudes of nodes 8 and 200 that were obtained for the following operating conditions: initial operating conditions (1, 5); operating conditions at reinforcement of weak ties (2, 6); operating conditions after the implementation of determined controls (3, 7); operating conditions under simultaneous implementation of the chosen controls and reinforcement of weak ties (4, 8).
An 11% decrease in the resistance of determined weak ties resulted in both a reduction in the mean square deviation of voltages at nodes 8 and 200 and a greater shift in the mean values with respect to the nominal voltage. Consequently, the probability that the voltage magnitudes at nodes 8 and 200 fall within the feasible limits decreased, as compared to the initial probability.
The implementation of control actions related to the change in the source reactive power at node 201 and to the regulation of transformation ratios of transformers, and the obtained shift in the mean values of voltage magnitudes at nodes 8 and 200 make it possible to increase the probability that they lie within the feasible limits.
Owing to the reduction in the mean square deviation and the shift in the mean, simultaneous implementation of control actions and reinforcement of weak ties increased the probability that the voltage magnitude at node 8 falls within the feasible limits. Graphs 2 and 4 in Figure 4, which correspond to the values of index ^AU and the probability that voltage magnitudes are within the
feasible limits, show that in this case both indices attest to the absence of critical voltage magnitudes.
Let us illustrate the operation of the algorithms with an example of the real electric network consisting of 207 nodes and 224 ties. The projection of the nodes and ties of this network on the plane in coordinates of the first and second singular vectors corresponding to nodal voltage magnitudes is presented in Figure 8. Voltage magnitudes at nodes 77, 78, 12, and 17 (110 kV) are sensor variables. The feasible interval for 110 kV voltages is taken equal to ±10 kV.
Figure 7. Probability density curves of changes in voltage magnitudes at nodes 8 and 200 for the initial operating conditions (1, 5), operating conditions after the reinforcement of weak ties (2, 6), operating conditions after the implementation of control actions (3, 7), and operating conditions under simultaneous implementation of control actions and reinforcement of weak ties (4, 8).
Figure 8. Projection of the real network graph on the plan in coordinates of the first (vi) and second (V2) right singular vectors that correspond to voltage magnitudes Table 2 shows that the same order of the nodes with sensor voltage magnitudes is obtained by the linear method of probabilistic load flow for the initial and end states. In the initial state the critical variable is only the voltage magnitude at node 12, whose probability of lying within the feasible limits is close to zero. After the shift of the mean, the end state with the 0.95 probability for
the critical value is determined. Such a result is obtained using only one control action found by the contribution factors method. The full control vector includes 93 components.
Table 2. Probabilistic characterictics of voltage magnitudes at the sensor nodes of real network for the initial (1) and end (2) states_
Nodes <* (kV) mv - Unom (kV) p
1 2 1 2 1 2
78 1.267 1.264 -5,69 -5,56 0,999 0,999
77 1,189 1,186 -4,24 -4,12 1 1
12 0,923 0,877 -11,96 -8,55 0,017 0,951
17 0,621 0,594 -7,21 -4,71 1 1
157 0,601 0,573 -7,41 -4,25 1 1
3 CONCLUSIONS
1. The methods of probabilistic load flow allow the detection of those sensor variables in the electric power system which can be detected on the basis of singular analysis.
2. A combination of the analytic probabilistic method with the scalar value of the first generalized disturbance is suggested to obtain probabilistic indices of variables in an inhomogeneous network.
3. An approach is proposed to solve the problem of selection of the control actions which provide the required probability that the controlled sensor variables fall within the feasible limits.
4. A method is suggested to search for a solution to the control problem with minimum number of controls on the basis of the data on tracing the power flows.
REFERENCES
[1] A.Z. Gamm, I.I. Golub, "Detection of weak places in a power system", Izv. RAS, Energetika, No.3, 1993, pp. 83-92, (in Russian).
[2] Gamm A.Z., Golub I.I., Observability of electric power systems, Novosibirsk: Nauka, 1990, p. 220, (in Russian).
[3] Fetzer E.E., Anderson P.M.,"Observability in the state estimation of power systems", IEEE Trans. on Power Apparatus and Systems, Vol. 94, No.6, 1975, pp. 1981-1988.
[4] Gamm A.Z., Golub I.I., "The problem of weak places in electric power systems", in Proc. 1995 International Conference on Electrical Power Engineering Power Tech, Stockholm, pp. 542-546.
[5] Gamm A.Z., Golub I.I., Bachry A., Styczynski Z.A.,"Solving several problems of power systems using spectral and singular analyses", IEEE Trans. Power Syst., Vol. 20, No.1, 2005, pp. 138-148.
[6] Bialek J., "Tracing the flow of electricity", IEE Generation Transmission &Distribution, Vol.143, No.4, 1996. pp. 313-320.
[7] Chen P., Chen Z., Bak-Jensen B., Probabilistic Load Flow: A Review http://vbn.aau.dk/ws/fbspretrieve/16272096/Probabilistic_load_flow_a_review.pdf
[8] Manusov V.Z., Mogirev V.V., Shepilov O.N., "Study of operating conditions of the USSR Unified Power System, considering the random character of initial data", Elektrichestvo, No. 10, 1983, pp. 3-6, (in Russian).
[9] Pugachev V.S., Theory of probabilities and mathematical statistics, M.: Nauka, 1979, p. 496, (in Russian).
[10] Borkowska B., "Probabilistic Load Flow", IEEE Trans. on Power Apparatus and Systems, Vol. 93, No. 3, 1974, pp. 752-755.
[11] . Kazakov I.E, Dostupov B.G., Statistical dynamics of nonlinear automatic systems, M.: Fizmatgiz, 1962, p. 332, (in Russian).
[12] Manusov V. Z., Lykin A.V., "Probabilistic analysis of steady states of electric power systems", Elektrichestvo, No.4, 1981, pp. 7-13, (in Russian).
[13] Chen X. Li, X., Yin X., Xiang T., H. Liu, The Algorithm of Probabilistic Load Flow Retaining Nonlinearity, Proc. of 2002 Power Con, Int. Conf. on Power System Technology, Kunming, Vol. 4, 2002, pp. 2111-2115.
[14] Su C. L., "Probabilistic load-flow computation using point estimate method", IEEE Trans. Power Systems, Vol. 20, No. 4, Nov. 2005, pp.1843-1851.
[15] Murashko N.A., Okhorzin Y.A., Krumm L.A. et al, Analysis and control of steady states of electric power systems, Novosibirsk: Nauka, 1987, p.240, (in Russian).
[16] Yudin D.B., Mathematical methods of control under the conditions of incomplete information, M.: Sov. radio, 1974. p.400, (in Russian).
[17] Gamm A.Z., Golub I.I., Application of Contribution factors at control of electric power system operating conditions, in Proc. 2005 International Conference on Electrical Power Engineering Power Tech, St. Petersburg, Rep. 649.
