A Stochastic Model for Determining Static Stability Margins in Electric Power Systems

Abstract

This paper aims to develop a method for determining margins of static aperiodic stability for electric power systems equipped with distributed generation plants. To this end, we used generalized equations of limiting modes in a stochastic formulation. Computer simulation showed that the developed methodology can be used in solving problems of operational control of the modes of electric power systems. On the basis of the results obtained, we arrived at the following conclusions: the modified equations do not allow the iterative process to converge to a trivial solution and, therefore, they ensure high reliability of results when determining stability margins in a stochastic statement; a technique based on the introduction of an additional variable can be used to improve the convergence of computational processes when determining the stability margins in a deterministic statement; the parameters of the limiting modes obtained in the deterministic and stochastic formulations may significantly differ; with an increase in the variance of the load graphs, the risk of stability violation significantly increases; at the same time, the amount of the margin determined on the basis of the Euclidean norm remains overly optimistic; in the illustrative example, a significant increase in the risk of stability violation takes place during planned and emergency shutdowns of the EPS elements.

1. Introduction

The design and operation of electric power systems (EPS) face problems of determining the limiting modes and margins of static aperiodic stability (SAS). Static stability is understood as the ability of the EPS to return to the steady state after small perturbations of the operating conditions, in which the changes in the parameters are very small compared to their average values. There are two types of static instability: aperiodic and dynamic. Static aperiodic stability is associated with a change in the active power balance in the EPS. Dynamic stability occurs when the automatic generator field regulators are set incorrectly. The first type of instability is associated with the appearance of positive real roots of the characteristic equation, and the second with the presence of complex roots with a positive real part.

The study in [1,2] develops the methods based on the equations of limiting modes proposed in [3,4]. In [5], the authors propose a non-iterative method to determine the SAS boundaries. The authors of [6] use power series to study the regions of existence of the EPS state. Reference [7] presents the results of the analysis of EPS states based on the tropical geometry of power balance equations over complex multipoles. The analysis of robust SAS is presented in [8]. Classical approaches to determining the static stability margins are given in [9,10]. Paper [11] presents the results of studies on stability margins in networks with DC lines. The research presented in [12] considers the stability margin problems for the networks with distributed generation. In [13], the authors focus on predicting the stability of electric power systems with a large share of wind generators.

In connection with the planned transition of the electric power industry to the technological platform of smart grids [14,15], which assumes a large-scale application of distributed generation (DG) plants [16,17], the formulated problems become relevant for power supply systems (PSS) since they determine the important factors associated with ensuring the required level of reliability [18].

Modern PSS tend to incorporate DG plants that use renewable energy sources, such as small-scale hydraulic power plants and wind farms that can be located at considerable distances from the consumption centers; at the same time, the margins of SAS are reducing, which has a significant effect on emergency control to ensure stability in post-emergency modes [18]. Therefore, the problem of building digital models for the operational determination of static stability margins in PSS with DG plants are of particular importance. One of the possible approaches to building such models is presented below. It is based on the equations of limiting modes [3] written in a stochastic formulation.

The analysis of SAS of electric power systems is carried out for the conditions of tuning automatic excitation regulators that eliminate the occurrence of self-swinging. The analysis procedure includes determination of the limiting modes, construction of the boundaries of the stability regions in the space of controlled parameters, and estimation of the SAS margin. The margin can be characterized by the roots of the characteristic equation based on finding the negative real root ${\alpha }_i$ located closest to the imaginary axis [10]. The margin can be estimated using the Jacobian of the steady-state mode equations that coincides with the free term of the characteristic polynomial. The values of synchronizing capacities $\frac{\partial P_i}{\partial {\delta }_i}$ can also be used as margin indicators.

It should be noted that most of the listed parameters are of little use for rationing of margins, since they characterize the PSS at the point of the mode under consideration. They are calculated using linearization of the system of equations, and the nonlinearity of power characteristics is not considered. Therefore, these indicators, borrowed from the theory of automatic control, do not find practical application due to the formality of estimates and insufficient clarity.

The articles published in 2020–2021 [19,20,21,22] consider the issues of static stability analysis and margin assessment in the framework of the intelligent control algorithms and active elements applied to increase SAS. For example, article [19] proposes improving the stability of multi-machine EPS by a static synchronous compensator controlled based on the ant colony algorithm. The study in [20] enhances the stability by creating an intelligent hybrid wind-solar farm as a static compensator. In [21], the problem of increasing SAS in a two-machine EPS is solved by using a fuel cell as a STATCOM. A real-time assessment of voltage stability in a large-scale EPS based on the spectrum estimation of the phasor measurement unit data is given in [22].

When determining the SAS margins, it is important to take into account the stochastic nature of changes in consumer loads, but the listed works do not address this issue.

The novelty of the results presented below is as follows:

• The stability margin is determined relying on the calculation of the limiting mode probability, which makes it possible to assess the risk of the planned measures to be taken to change the mode.
• The essential difference between the approach proposed below and the known methods for determining limiting modes is that the corresponding Jacobian matrix is nondegenerate at the solution point, which ensures reliable convergence of computational processes when calculating stability margins [3].
• Reliable result is provided and the iterative process does not converge to the point of a trivial solution corresponding to the initial steady state.

Of particular relevance is the development of a method to stochastically determine stability margins in the isolated power systems [23] involving renewable energy sources. Such systems can be established by applying a cyber-physical approach [24,25] based on the most advanced information and computer technologies, such as artificial intelligence, big data, the Internet of things, quantum computing, and others. However, the core of the virtual part of the cyber-physical power system should rely on digital models based on improved algorithms designed to solve traditional electric power problems, which include the determination of SAS margins as well.

2. Problem Statement

In the problems of operational control and determination of reliability of the EPS, as well as when choosing the control actions of emergency automation, high-speed high-validity and reliability methods for assessing margins are required. The validity and reliability requirements are met by the SAS margin indicator that corresponds to the distance in the space of controlled parameters Y from the point Y₀ of the analyzed mode to the limiting hypersurface $L_W$. Assessment of the SAS margin can be carried out using the vector

$$K={\left[k_1{k}_2\dots k_m\right]}^T,$$

Which contains margin coefficients with respect to parameters $y_i$: noindent.

$$k_i={\mathsf{\mu }}_i\left(y_{iL}-y_{i0}\right);{\mu }_i=\frac{1}{k_{inorm}{y}_{i0}},$$

where $y_{i0}, y_{iL}$ are the values of the i-th parameter corresponding to the initial and limiting modes; $k_{inorm}$ are normative coefficients.

The critical (most dangerous) direction of loading can be found by solving the following optimization problems:

$$\min Z=\underset{X_L}{\min } \underset{i}{\max }\left|k_i\left(X_L\right)\right|$$

(1)

or

$$\mathrm{minZ}=\underset{X_L}{\min }\sqrt{{\displaystyle \sum _{i\ne j}{k}_i^2\left(X_L\right)}},$$

(2)

where $X_L$ is a vector of dependent variables at the points of the limit surface; $Z$ is the margin value.

Condition (1) does not fully characterize the closeness of the mode to the boundary $L_w$, because it does not take into account changes in all variables $y_j$ included in the loading vector. The criterion defined by the Euclidean norm (2) does not have this shortcoming and corresponds to the shortest distance from point Y₀ to $L_w$ in space K.

Further development of methods for determining the SAS margins based on the search for the critical direction of load increase can be carried out using the stochastic approach described below.

3. Stochastic Approach to Determining the Stability Margin

The algorithm that implements the stochastic approach to determine the stability margins is based on a generalization of the equations of limiting modes [3].

Components of the vector $Y$ that represent active and reactive powers of generators and loads are random variables. With irregular changes in these parameters, it is possible to reach the boundary $L_w$; in this case, reliable operation of the PSS will be ensured if the scattering hyperellipsoid (Figure 1) [4]

$$D{Y}^T{S}^{-1}DY-C_L^2=0$$

(3)

With the center at the point $MY$ does not have common points with the hypersurface $L_W$. Here $DY=Y-MY$; $M$ denotes mathematical expectation; $S=M\left[\left(Y-MY\right){\left(Y-MY\right)}^T\right]$ is a covariance matrix; $C_L$ is a value that determines the probability of finding the parameters of Y inside the hyperellipsoid (3). The covariance matrix S participating in (3) can be found using known methods for assessing the state of the PSS [26].

The stochastic approach to assess the margin of SAS can be formulated as follows:

Find

$$C_L=\min {\left(D{Y}^T{S}^{-1}DY\right)}^{\frac{1}{2}}$$

With constraints in the form of the steady-state equations

$$F\left(X,MY+DY\right)=0,$$

where F is a nonlinear vector function; X is unregulated operating parameters, which are the magnitudes and phases of nodal voltages or their real and imaginary components.

To solve the problem of minimizing $C_L$, build a Lagrange function

$$L\left(X,MY+DY,\Lambda \right)={\left(D{Y}^T{S}^{-1}DY\right)}^{\frac{1}{2}}+F^T\left(X,MY+DY\right)\Lambda ,$$

where $\Lambda$ is a vector of undefined factors.

The following conditions correspond to the minimum of L:

$$\begin{array}{l}\frac{\partial L}{\partial DY}=S^{-1}DY{\left(D{Y}^T{S}^{-1}DY\right)}^{-\frac{1}{2}}+{\left(\frac{\partial F}{\partial DY}\right)}^T\Lambda =0;\\ \frac{\partial L}{\partial X}={\left(\frac{\partial F}{\partial X}\right)}^T\Lambda =0;\\ \frac{\partial L}{\partial \Lambda }=F\left[X,MY+DY\right]=0,\end{array}\},$$

(4)

where $\frac{\partial F}{\partial X}$ is a Jacobian of the steady-state mode equations; $\frac{\partial F}{\partial DY}$ is a block diagonal matrix.

It should be taken into consideration that Equation (4) has two solutions:

• The trivial solution corresponds to zero values of the vectors

  $$\Lambda ={\left[0\dots 0\right]}^T, DY={\left[0\dots 0\right]}^T.$$
• The desired solution exists if one or several components of the vectors $\Lambda$ and $DY$ are not equal to zero; at the same time, the equations $\left(\frac{\partial F}{\partial X}\right){\Lambda }^T=0$ correspond to the condition

  $$\det \left(\frac{\partial F}{\partial X}\right)=0,$$

Hence, the nontrivial solution corresponds to the hypersurface of the limiting modes.

Equation (4) can written as

$$\begin{array}{l}\frac{\partial L}{\partial DY}=S^{-1}DY+{\left(\frac{\partial F}{\partial DY}\right)}^T\Lambda {\left(D{Y}^T{S}^{-1}DY\right)}^{\frac{1}{2}}=0;\\ \frac{\partial L}{\partial X}=\left(\frac{\partial F}{\partial X}\right)\Lambda =0;\\ \frac{\partial L}{\partial \Lambda }=F\left[X,MY+DY\right]=0.\end{array}\}$$

Introduce a change of variables

$$R={\left(D{Y}^T{S}^{-1}DY\right)}^{\frac{1}{2}}\Lambda ,$$

Then, system (4) can be written in the following form

$$\begin{array}{l}{S}^{-1}DY+{\left(\frac{\partial F}{\partial DY}\right)}^{T}R=0;\\ {\left(\frac{\partial F}{\partial X}\right)}^{T}R=0;\\ \\ F\left(X,MY+DY\right)=0.\end{array}\},$$

(5)

From system (5), we can find a vector

$$DY=-S{\left(\frac{\partial F}{\partial DY}\right)}^{T}R=0$$

Substituting the resulting value of $DY$ into the third equation, we can transform the system to the form:

$$\begin{array}{l}F\left[X,MY-S{\left(\frac{\partial F}{\partial DY}\right)}^{T}R\right]=0;\\ {\left(\frac{\partial F}{\partial X}\right)}^{T}R=0.\end{array}\}$$

(6)

The results of computer simulation show that the iterative process for Equation (6) may converge to a trivial solution. To eliminate this, we can employ the following technique. In system (6), we introduce an additional equation, which corresponds to a nonzero length of the vector $R$, and a variable γ, which ensures the “balancing” of the first vector equation:

$$\begin{array}{l}F\left(X,R,\gamma \right)=F\left[X,MY-\gamma S{\left(\frac{\partial F}{\partial DY}\right)}^{T}R\right]=0;\\ V\left(X,R\right)={\left(\frac{\partial F}{\partial X}\right)}^{T}R=0;\\ \\ w\left(R\right)=R^{T}R-1=0.\end{array}\}$$

(7)

A similar technique can be used to improve the efficiency of determining the stability margin in a deterministic statement [3].

Equation (7) is nonlinear and can be solved only by iterative methods, for example, Newton’s method. In this case, the following system of linear equations is solved at each interaction:

$$\left[\begin{array}{ccc}\frac{\partial F}{\partial X}& \gamma S{\left(\frac{\partial F}{\partial DY}\right)}^T& S{\left(\frac{\partial F}{\partial DY}\right)}^{T}R\\ \frac{\partial V}{\partial X}& {\left(\frac{\partial F}{\partial X}\right)}^T& 0\\ 0& 2R^T& 0\end{array}\right]\left[\begin{array}{c}\begin{array}{l}\Delta X\\ \end{array}\\ \begin{array}{l}\Delta R\\ \end{array}\\ \Delta \gamma \end{array}\right]=-\left[\begin{array}{c}\begin{array}{l}F\\ \end{array}\\ \begin{array}{l}V\\ \end{array}\\ w\end{array}\right]$$

Elements of matrix $\frac{\partial V}{\partial X}$ included in this system are determined by the expressions

$$\frac{\partial \mathbf{v}}{\partial \mathbf{X}}=\frac{\partial }{\partial \mathbf{X}}\left[{\left(\frac{\partial \mathbf{F}}{\partial \mathbf{X}}\right)}^T\mathbf{R}\right]=\frac{\partial }{\partial \mathbf{X}}{\left[\left({\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_1}{r}_i}\right)\left({\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_1}{r}_i}\right)\dots \left({\displaystyle \sum _{i=1}^n\frac{\partial f_i}{\partial x_1}{r}_i}\right)\right]}^T$$

By differentiating, one can obtain

$$\frac{\partial \mathbf{V}}{\partial \mathbf{X}}=\left[\begin{array}{cccc}{\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_1\partial x_1}{r}_i}& {\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_1\partial x_2}{r}_i}& \dots & {\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_1\partial x_n}{r}_i}\\ {\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_2\partial x_1}{r}_i}& {\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_2\partial x_2}{r}_i}& \dots & {\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_2\partial x_n}{r}_i}\\ \dots & \dots & \dots & \dots \\ {\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_n\partial x_1}{r}_i}& {\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_n\partial x_2}{r}_i}& \dots & {\displaystyle \sum _{i=1}^n\frac{{\partial }^2{f}_i}{\partial x_n\partial x_n}ri}\end{array}\right]$$

Hence

$$\frac{\partial \mathbf{V}}{\partial \mathbf{X}}={\displaystyle \sum _{i=1}^n{\Gamma }_i{r}_i},$$

where ${\Gamma }_i$ is Hessian matrix of function f₁.

If the steady-state equations can have the form

$$F\left(X,Y\right)=Y-F\left(X\right)=0,$$

then the following equality is true

$${\left(\frac{\partial F}{\partial DY}\right)}^T=E$$

where E is an identity matrix. Such a situation occurs when loads are specified by constant power take-offs.

With an implicit Y dependency of X and by using the Cartesian coordinate system, one can write

$$f_{2i-1}\left(X,Y\right)=P_{i0}+d{P}_i\left(U_i^{\prime },U_i^{″},\right)-P_{ci}\left(\begin{array}{ccccc}{U}_1^{\prime }& U_1^{″}& \dots & U_n^{\prime }& U_n^{″}\end{array}\right)=0;$$

$$f_{2i}\left(X,Y\right)=Q_{i0}+d{Q}_i\left(U_i^{\prime },U_i^{″}\right)-Q_{ci}\left(\begin{array}{ccccc}{U}_1^{\prime }& U_1^{″}& \dots & U_n^{\prime }& U_n^{″}\end{array}\right)=0;$$

Or

$$f_{2i}\left(X,Y\right)=U_{Zi}+D{U}_{Zi}-\sqrt{{\left(U_i^{\prime }\right)}^2+{\left(U_i^{″}\right)}^2=0};$$

where $P_{i0}=P_{Gi0}-P_{Hi0}$, $Q_{i0}=Q_{Gi0}-Q_{Hi0}$ are active and reactive power injections at node i in the initial mode; $P_{Gi0},Q_{Gi0}$ are active and reactive powers of generators; $P_{Hi0},Q_{Hi0}$ are active and reactive powers of loads; $U_i^{\prime },U_i^{″}$ are real and imaginary components of nodal voltages; $U_{\mathrm{Z}i}$ is a set magnitude of voltage at the i-th node; $d{P}_i$, $d{Q}_i$, $D{U}_{Zi}$ are components of vector $DY$; $P_{Ci},Q_{Ci}$ are power flows from node i to the network.

To obtain $\frac{\partial F}{\partial DY}$ at an implicit X dependency of Y, steady-state conditions should be written in an expanded form:

For generator nodes

$$f_{2i-1}\left(X,Y\right)=P_{Gi}+D{P}_{Gi}-\left(P_{Hi}+D{P}_{Hi}\right)h_{1i}-P_{Ci}\left(U_1^{\prime },U_1^{″},\dots U_n^{\prime },U_n^{″}\right)=0;$$

(8)

$$f_{2i}\left(X,Y\right)=U_{Zi}+D{U}_{Zi}-\sqrt{{\left(U_i^{\prime }\right)}^2+{\left(U_i^{″}\right)}^2=0};$$

For load nodes

$$f_{2i-1}\left(X,Y\right)=P_{Gi}+D{P}_{Gi}-\left(P_{Hi}+D{P}_{Hi}\right)h_{1i}-P_{Ci}\left(U_1^{\prime },U_1^{″},\dots U_n^{\prime },U_n^{″}\right)=0;$$

(9)

$$f_{2i}\left(X,Y\right)=-\left(Q_{Hi}+D{Q}_{Hi}\right)h_{2i}-Q_{Ci}\left(U_1^{\prime },U_1^{″},\dots U_n^{\prime },U_n^{″}\right)=0,$$

where $D{P}_{Gi}$, $D{U}_{Zi}$, $D{P}_{Hi}$, $D{Q}_{Hi}$ are components of vector $DY$;

$$h_{1i}=\left[a_{0i}+a_{1i}\left(\frac{U_i}{U_{NOMi}}\right)+a_{2i}{\left(\frac{U_i}{U_{NOMi}}\right)}^2\right];h_{2i}=\left[b_{0i}+b_{1i}\left(\frac{U_i}{U_{NOMi}}\right)+b_{2i}{\left(\frac{U_i}{U_{NOMi}}\right)}^2\right];$$

$a_{0i};a_{1i};a_{2i};b_{0i};b_{1i};b_{2i}$ are coefficients of polynomials approximating the static characteristics of the load; $U_{NOMi}$ is rated voltage of the i-th node.

For the nodes corresponding to Equations (8) and (9), matrix $\frac{\partial F}{\partial DY}$ will include the following parts:

For generator nodes

Components of Vector DY	$f_{2i-1}\left(X,Y\right)$	$f_{2i}\left(X,Y\right)$
$D{P}_{Gi}$	$\frac{\partial f_{2i-1}}{\partial D{P}_{Gi}}$	0
$D{P}_{Hi}$	$\frac{\partial f_{2i-1}}{\partial D{P}_{Hi}}$	0
$D{U}_{Zi}$	0	$\frac{\partial f_{2i}}{\partial D{U}_{Zi}}$

For load nodes

Components of Vector DY	$f_{2i-1}\left(X,Y\right)$	$f_{2i}\left(X,Y\right)$
$D{P}_{Hi}$	$\frac{\partial f_{2i-1}}{\partial D{P}_{Hi}}$	0
$D{Q}_{Hi}$	0	$\frac{\partial f_{2i}}{\partial D{Q}_{Hi}}$

where $\frac{\partial f_{2i}}{\partial D{U}_{Zi}}$ = $\frac{\partial f_{2i-1}}{\partial D{P}_{Gi}}$ = 1; $\frac{\partial f_{2i-1}}{\partial D{P}_{Hi}}$ = $h_{1i}$; $\frac{\partial f_{2i}}{\partial D{Q}_{Hi}}$ = $h_{2i}$.

Based on the value $C_L$ obtained by solving Equation (7), one can calculate $R_{SAS}$ that determines probabilistic assessment of the risk of SAS loss; in this case, an approach similar to that proposed in [27] can be used to assess financial risks

$$R_{SAS}=\exp \left(-\frac{D{Y}_L^T{S}^{-1}D{Y}_L^{}}{2}\right)=\exp \left(-\frac{C_L^2}{2}\right)$$

where $D{Y}_L^{}$ is the value of vector $DY$ at the point of contact between the hyperellipsoid (3) and the hypersurface L_W.

System (7), compared to the steady-state equations that form the first group of equations of this system, has an increased dimension. However, the exclusion of multistep procedures of discrete weighting [10] ensures high efficiency of calculations and reliability of results in the problems of operational control of electric power systems.

4. Simulation Results

As an example, below are the results of determining the stability margins for a three-node EPS model shown in Figure 2. The model corresponds to a 110 kV network supplying two load nodes equipped with distributed generation (DG) plants with a rated power of 50 MW. Automatic regulators of DG plants provide voltage stabilization on 110 kV busbars. Figure 3 shows the graphs of the power injection changes

$$P_k\left(t\right)=P_{Hk}\left(t\right)-P_{Gk},$$

where $P_{Hk}\left(t\right)$ are consumer loads with a random time law; $P_{GK}$—capacities of generators of the DG plants, which are assumed to be constant; k = 1, 2.

The approach used in this example was based on the following assumptions:

• a stochastic approach was used for consumer loads; at the same time, their capacities changed according to the graphs shown in Figure 2;
• deterministic models traditionally used in practice were employed for the powers of the generator; it should be noted that the algorithm described above can also be used without modifications in the case of a stochastic nature of changes in the capacities of generators.

Table 1. Dependence of the parameters of the limiting mode on the standard deviation $\sigma$ of the load capacities.

$\mathit{\sigma },\%$	${\mathit{C}}_{\mathit{L}}, \mathbf{p}.\mathbf{u}.$	${\mathit{R}}_{\mathit{S}\mathit{A}\mathit{S}}, \%$	${\mathit{y}}_{\mathit{L}1}^{\left(\mathit{S}\right)}, \mathbf{MW}$	${\mathit{y}}_{\mathit{L}2}^{\left(\mathit{S}\right)}, \mathbf{MW}$	Z, %
10	6.4	1.28·10–7	–248.00	–275	74.76
20	4.03	0.03	–273.00	–245	63.67
30	2.5	4.4	–276.00	–241	62.55
40	1.9	16.4	–278.00	–239	62.19

Note: the values of Z were calculated using formula (2).

presents simulation results obtained on the basis of an experimental code implemented in Mathcad. Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 provide the corresponding illustrations. Table 1 shows the values of the parameters $y_{L1}^{\left(S\right)}$,$y_{L2}^{\left(S\right)}$ characterizing the limiting modes obtained at different values of the standard deviation $\sigma$ of the consumer load graphs. In addition, Table 1 provides values of $C_L$ corresponding to hyperellipsoid (3), and the values that determine the probabilistic risk assessments of the SAS loss.

It follows from Table 1 that, with significant load variations, the probability of stability loss can be quite high and reach 16%. However, the stability margin determined by the Euclidean norm (2), which does not take into account the stochastic nature of the operating conditions, will be even higher (62%). Thus, the margin value determined on the basis of the Euclidean norm may be overly optimistic.

Figure 4 shows the results of determining the limiting modes in a stochastic statement. For comparison, we performed calculation of a limiting mode using a deterministic approach; it was assumed that

$$S=E,$$

where $E$ is an identity matrix. The results of calculating the limiting mode in a deterministic statement are illustrated in more detail in Figure 5. As seen in Figure 5, the limiting mode parameters obtained in deterministic $Y_L^{\left(D\right)}$ and stochastic $Y_L^{\left(S\right)}$ statements can differ greatly. Figure 6 shows the results of determining the limiting mode for σ = 40%.

Figure 7 shows the results of determining the limiting mode when the 50 MW generator of the DG plant is disabled at node 1. Figure 8 provides illustration for the results obtained when one circuit of the electric transmission line connecting nodes 2 and 3 is disabled.

Using Table 1, we drew graphs for $C_L^{}=C_L^{}\left(\sigma \right)$, $R_{SAS}=R_{SAS}\left(\sigma \right)$, and $Z=Z\left(\sigma \right)$ shown in Figure 9 and Figure 10.

It should be noted that at the solution to Equation (7), vectors of the normals to the limit hypersurface $L_w$, subject to

$$\det \left(\frac{\partial F}{\partial X}\right)=0,$$

and to hyperellipsoid (3) are collinear.

It was shown in [5] that the vector $R$ determined by solving Equation (7) coincides with the direction of the normal to $L_w$. The direction of the normal to the hyperellipsoid (3) can be found by calculating the gradient vector

$$N=grad \left(D{Y}^T{S}^{-1}DY-C_L^2\right)=\left({\left(S^{-1}\right)}^T+S^{-1}\right)DY$$

Collinearity of vectors $R$ and $N$ can be determined by finding the angle between them:

$$\cos \left(N,R\right)=\frac{N^{T}R}{\sqrt{\left(N^{T}N\right)}\cdot \sqrt{\left(R^{T}R\right)}}$$

Calculations of the limiting mode obtained for $\mathsf{\sigma }$ = 15% provide the following results:

$$y_{L1}^{\left(S\right)}=-261 \mathrm{MW}; y_{L2}^{\left(S\right)}=-260 \mathrm{MW}; R={\left(\begin{array}{cc}0.767& 0.641\end{array}\right)}^T;$$

$$N={\left(\begin{array}{cc}-0.349& -0.291\end{array}\right)}^T; \cos \left(N,R\right)=-1.$$

These results confirm the theoretical findings about collinearity of the vectors R and N.

5. Conclusions

Based on the results of computer simulation, the following conclusions can be made:

• The modified equations that do not allow the iterative process to converge to the trivial solution ensure high reliability of results when determining the stability margins in a stochastic statement; a technique based on the introduction of an additional variable can be used to improve the convergence of computational processes in determining the stability margins in a deterministic statement.
• The parameters of the limiting modes obtained in the deterministic and stochastic formulations may significantly differ.
• With an increase in the variance of load graphs, the risk of stability violation increases significantly; at the same time, the amount of the margin determined on the basis of the Euclidean norm remains overly optimistic.
• In the illustrative example, a significant increase in the risk of stability violation occurs during planned and emergency shutdowns of the EPS components.
• At present, the concept of cyber-physical systems is being actively developed. It provides the foundation for building modern, efficient, and reliable EPSs. This concept presupposes the use of not only physical objects, but also deep integration of digital components that process a large amount of information coming from sensors, which is then used to simulate and control physical objects. The proposed method for determining the stability margins in a stochastic statement can be effectively used to solve the problems of analyzing SAS in cyber-physical EPS equipped with distributed generation units.
• Further research could be aimed at factoring in the asymmetry of loads when determining the stability margins in a stochastic formulation based on the use of phase coordinates to write the equations of limiting modes [28].