TOWARD A UNIFIED ALGORITHM FOR CALCULATING STATIC STABILITY AND STEADY-STATE CHANGES IN AN ELECTRIC POWER SYSTEM

R. I. BORISOV, V. I. GOTMAN, Yu. V. KHRUSHCHEV

(Presented by the Department of Electrical Systems and Networks)

When a stationary disturbance is present in a system, it becomes necessary to solve a number of problems.

1. Determination of the parameters of the new steady state (calculation of the change in operating conditions).

2. Calculation of static stability:

a) verification of a stable transition from one steady state to another;

b) determination of the parameter values of the limiting state in order to assess the static-stability margin of the new state.

Usually the indicated problems are considered as separate, independent ones and are described by different equations.

The operating state of an electric power system is described by equations of nodal voltages or loop currents, in which the loads are represented by constant resistances [1]. This method of accounting for the load is only a certain approximation of its operating properties. A change in the magnitude of the load leads to the need to recalculate the elements of the matrix of self and mutual admittances, or to correct them, which entails certain inconveniences [2]. In static-stability calculations, equations of nodal powers are used, with the same assumptions regarding the accounting of the load [3].

However, the problems of calculating the change in the operating state of an electric power system and determining the limiting values of power flows over sections under static-stability conditions can be described by a unified algorithm; in this case loads and generators are taken into account by their static characteristics.

The equations given below make it possible to solve the following problems:

1. To determine the operating-state parameters of an electric power system in the presence of control and disturbing actions specified in the form of changes in transformation ratios, active and reactive powers of loads and generators, voltage or emf, as well as in the presence of symmetrical (three-phase) network disturbances, and to check the static stability of the new state according to practical criteria.

2. To calculate the coefficients of transfer functions for any operating-state parameters at constant frequency.

3. To equivalent sections of a complex electric power system by generalized static characteristics and to find the coefficients

the curvature of the equivalent part when it is represented by a two-terminal network.

1. Calculate the limiting operating conditions of the system corresponding to practical stability criteria.

For an electric system, two groups of equations are compiled: the first group describes the operating condition of the electric network, represented by a passive multi-terminal network with linear values of self and mutual conductances; the second describes the static characteristics of the generating sources and loads connected to the corresponding nodes (vertices) of the multi-terminal network (Fig. 1).

The desired increment equations are obtained by expanding the original equations in a Taylor series (retaining only the linear terms of the expansion) and by their subsequent joint solution. The equations are solved with respect to specified unknowns for an arbitrary combination of the above-listed disturbances; at the same time, any constraints on operating parameters may be taken into account.

Fig. 1. Diagram of the electric system

Calculation of Changes in the Steady-State Operating Condition of the System

The vertices of the passive multi-terminal network are taken to be all nodes of the electric network to which generators, synchronous compensators, loads, and branches of transverse conductance are connected. A regulating transformer is introduced into the equivalent circuit in the form of an independent branch, the ends of which are also considered separate nodes. If none of the above-mentioned elements is connected to some node, then the corresponding vertex of the multi-terminal network is regarded as the point of application of a fictitious (zero) load.

To calculate increments of the operating parameters, the following must be specified: the circuit and parameters of the system, the initial operating condition (voltage vectors at the nodes), and the equations of the static characteristics of the loads or their voltage-slope coefficients.

The operating condition of a complete multi-terminal network containing \(n\) nodes and \(m\) branches can be described by \(2n\) equations of nodal powers. These equations, for each node \(i\), have the form

\[ P_i=U_i^2 g_{ii}+U_i\sum_{\substack{j=1\\ j\ne i}}^n U_j y_{ij}\sin(\delta_i-\delta_j-\alpha_{ij}). \tag{1} \]

\[ Q_i=U_i^2 b_{ii}-U_i\sum_{\substack{j=1\\ j\ne i}}^n U_j y_{ij}\cos(\delta_i-\delta_j-\alpha_{ij}). \tag{2} \]

The positive directions of power are taken toward the nodes. The powers of sources are taken with a plus sign, and those of receivers with a minus sign.

In equations (1), (2), the voltages of nodes adjacent to node \(i\), and the admittances of the branches incident to it, \(y_{ij}=q_{ij}-j b_{ij}\), are referred to the transformation ratio of node \(i\) (see the Appendix). Thus the transformation coefficients explicitly enter the equations. The admittances of unregulated transformers and of branches incident to them may be referred to some base voltage.

Equations (1), (2) are functions of the magnitudes and phases of the node voltages and of the transformation coefficients of regulated transformers. Expanding the indicated equations in a Taylor series with respect to the independent variables \(U,\delta,K_{\mathrm{T}}\), we obtain the increments of nodal powers, which in matrix form are written as

\[ \left\| \begin{array}{c} \Delta P\\ \Delta Q \end{array} \right\| = \left\| \begin{array}{ccc} \dfrac{\partial P}{\partial U} & \dfrac{\partial P}{\partial \delta} & \dfrac{\partial P}{\partial K_{\mathrm{T}}}\\[6pt] \dfrac{\partial Q}{\partial U} & \dfrac{\partial Q}{\partial \delta} & \dfrac{\partial Q}{\partial K_{\mathrm{T}}} \end{array} \right\| \left\| \begin{array}{c} \Delta U\\ \Delta \delta\\ \Delta K_{\mathrm{T}} \end{array} \right\|. \tag{3} \]

The square submatrices \(\partial P/\partial U\) and \(\partial P/\partial \delta\) have the form

\[ \frac{\partial P}{\partial U} = \left\| \begin{array}{cccc} \dfrac{\partial P_1}{\partial U_1} & \dfrac{\partial P_1}{\partial U_2} & \cdots & \dfrac{\partial P_1}{\partial U_n}\\ \cdot & \cdot & \cdots & \cdot\\ \dfrac{\partial P_n}{\partial U_1} & \dfrac{\partial P_n}{\partial U_2} & \cdots & \dfrac{\partial P_n}{\partial U_n} \end{array} \right\|, \qquad \frac{\partial P}{\partial \delta} = \left\| \begin{array}{cccc} \dfrac{\partial P_1}{\partial \delta_1} & \dfrac{\partial P_1}{\partial \delta_2} & \cdots & \dfrac{\partial P_1}{\partial \delta_n}\\ \cdot & \cdot & \cdots & \cdot\\ \dfrac{\partial P_n}{\partial \delta_1} & \dfrac{\partial P_n}{\partial \delta_2} & \cdots & \dfrac{\partial P_n}{\partial \delta_n} \end{array} \right\|. \tag{4} \]

The submatrices \(\partial Q/\partial U\), \(\partial Q/\partial \delta\) have an analogous form. The rectangular submatrix \(\partial P/\partial K_{\mathrm{T}}\) is analogous to \(\partial Q/\partial K_{\mathrm{T}}\) and is written as

\[ \frac{\partial P}{\partial K_{\mathrm{T}}} = \left\| \begin{array}{cccc} \dfrac{\partial P_1}{\partial K_{\mathrm{T}1}} & \dfrac{\partial P_1}{\partial K_{\mathrm{T}2}} & \cdots & \dfrac{\partial P_1}{\partial K_{\mathrm{T}m}}\\ \cdot & \cdot & \cdots & \cdot\\ \dfrac{\partial P_n}{\partial K_{\mathrm{T}1}} & \dfrac{\partial P_n}{\partial K_{\mathrm{T}2}} & \cdots & \dfrac{\partial P_n}{\partial K_{\mathrm{T}m}} \end{array} \right\|. \tag{5} \]

In each row of the indicated submatrices, the nonzero elements are the derivatives with respect to the transformation coefficients of transformers with OLTCs included in branches incident to node \(i\). The partial derivatives found from (1) and (2) are given in the Appendix.

The column submatrices \(\Delta U, \Delta\delta, \Delta P, \Delta Q\) of order \(n\) contain the increments of the modulus, phase, and power at all nodes of the network. The column submatrix \(\Delta K_{\mathrm{T}}\) of order \(m\) contains the increments of the transformation ratios for all branches; moreover, for branches in which there is no transformer, or \(K_{\mathrm{T}}=\mathrm{const}\), \(\Delta K_{\mathrm{T}}=0\).

For generality of presentation, we assume that the transformers provide longitudinal-transverse voltage regulation and, consequently, the transformation ratios are complex quantities \([1, 4]\). The increment of the modulus of the transformation ratio in (3) can be expressed in terms of increments of the components of the vector \((K'_{\mathrm{T}}\) and \(K''_{\mathrm{T}})\) from the equation

\[ K_{\mathrm{T}i}=\sqrt{K_{\mathrm{T}i}^{\prime 2}+K_{\mathrm{T}i}^{\prime\prime 2}}, \tag{6} \]

whence

\[ \Delta K_{\mathrm{T}i}=\mathring K'_{\mathrm{T}i}\Delta K'_{\mathrm{T}i}+\mathring K''_{\mathrm{T}i}\Delta K''_{\mathrm{T}i}, \tag{7} \]

where

\[ \mathring K'_{\mathrm{T}i}=K'_{\mathrm{T}i}/K_{\mathrm{T}i}, \qquad \mathring K''_{\mathrm{T}i}=K''_{\mathrm{T}i}/K_{\mathrm{T}i}. \tag{8} \]

The matrix \(\Delta K_{\mathrm{T}}\) is equal to

\[ \Delta K_{\mathrm{T}}=\left\|\mathring K'_{\mathrm{T}}\ \mathring K''_{\mathrm{T}}\right\| \begin{Vmatrix} \Delta K'_{\mathrm{T}}\\ \Delta K''_{\mathrm{T}} \end{Vmatrix}, \tag{9} \]

where the submatrices \(\mathring K'_{\mathrm{T}}, \mathring K''_{\mathrm{T}}\) are diagonal.

The elements connected to the nodes of the multiport network (loads, generators, synchronous compensators) may be described by their static characteristics.

The equation of the static load characteristic with respect to voltage has the form \([7]\)

\[ S_i(U_i)=P_i(U_i)+jQ_i(U_i) \tag{10} \]

and, correspondingly, the slope coefficients are

\[ \alpha_{i\mathrm{n}}=dP_{i\mathrm{n}}/dU_i, \qquad \beta_{i\mathrm{n}}=dQ_{i\mathrm{n}}/dU_i. \tag{11} \]

The static characteristics of generators may be represented in the form \([5]\)

\[ Q_{i\mathrm{g}}=Q(U_i,P_{i\mathrm{g}}). \tag{12} \]

The slope coefficients are determined with allowance for automatic excitation regulators at fixed settings of the amplifying elements,

\[ \beta_{i\mathrm{g}}=\partial Q_{i\mathrm{g}}/\partial U_i\,(P_{i\mathrm{g}}=\mathrm{const}), \qquad \eta_{i\mathrm{g}}=\partial Q_{i\mathrm{g}}/\partial P_{i\mathrm{g}}\,(U_i=\mathrm{const}). \tag{13} \]

The static characteristics of synchronous compensators are

\[ Q_{i\mathrm{k}}=Q(U_i) \tag{14} \]

and the slope coefficient, under the same conditions as for generators, is

\[ \beta_{i\mathrm{k}}=dQ_{i\mathrm{k}}/dU_i. \tag{15} \]

In what follows we assume that the values of all coefficients \(\alpha, \beta, \eta\) in the vicinity of the initial operating condition are known.

The power increments at the nodes of the network in equation (3) can be expressed through the power increments of generators, loads, and synchronous compensators. The total power increments of each two-terminal element \(i\) can be represented in the form of specified increments \((\Delta P_{i\mathrm{z}}, \Delta Q_{i\mathrm{z}})\), which, together with \(\Delta K_{\mathrm{T}}\), constitute independent variables—control and disturbance actions and increments that—

which are forced in nature and are a function of the above-mentioned influences, varying in accordance with their static characteristics. These forced components can be expressed through \(\Delta U\), \(\Delta P_{\mathrm{g}}\) and the slope coefficients of the static characteristics, \(\alpha, \beta, \eta\).

Proceeding from the foregoing, the matrix of total power increments of the two-terminal networks, taking (10)--(15) into account, can be written as

\[ \left\|\begin{array}{c} \Delta P\\ \Delta Q \end{array}\right\|_{\mathrm{l}} = \left\|\begin{array}{c} \Delta P_{\mathrm{s}}\\ \Delta Q_{\mathrm{s}} \end{array}\right\| + \left\|\begin{array}{c} \Delta P_{b}\\ \Delta Q_{b} \end{array}\right\| = \left\|\begin{array}{c} \Delta P_{\mathrm{s}}\\ \Delta Q_{\mathrm{s}} \end{array}\right\| + \left\|\begin{array}{cc} \alpha & 0\\ \beta & \eta \end{array}\right\| \left\|\begin{array}{c} \Delta U\\ \Delta P_{\mathrm{s}} \end{array}\right\|. \tag{16} \]

Here \(\Delta P_{\mathrm{s}}, \Delta Q_{\mathrm{s}}\) are column submatrices of the specified disturbances; \(\alpha, \beta, \eta\) are diagonal submatrices of the slope coefficients of the static characteristics, of order \(n\).

Replacing the power increments at the network nodes in (3) by the equal power increments of the two-terminal networks from (16), after transformations taking (9) into account, we obtain the desired matrix equation for \(2n\) increments (\(n\) increments \(\Delta U\), \(n - 1\) increments \(\Delta \delta\), and \(\Delta P_n\)) for given \(2n + m - 1\) disturbances \(\Delta P_{\mathrm{s}}, \Delta Q_{\mathrm{s}}, \Delta K_{\mathrm{T}}\) (including disturbances equal to zero):*

\[ \left\|\begin{array}{cc} W_{pU} & W_{p\delta}\\ W_{qU} & W_{q\delta} \end{array}\right\| \left\|\begin{array}{c} \Delta U\\ \Delta\delta_{(n)} \end{array}\right\| = \left\|\begin{array}{c} G_p\\ G_q \end{array}\right\|. \tag{17} \]

The column submatrices of disturbances \(G_p, G_q\) have the form

\[ G_p = \Delta P_{\mathrm{s}} - \frac{\partial P}{\partial K_{\mathrm{T}}'}\Delta K_{\mathrm{T}}' - \frac{\partial P}{\partial K_{\mathrm{T}}''}\Delta K_{\mathrm{T}}'', \tag{18} \]

\[ G_q = \Delta Q_{\mathrm{s}} + \eta \Delta P_{\mathrm{s}} - \frac{\partial Q}{\partial K_{\mathrm{T}}'}\Delta K_{\mathrm{T}}' - \frac{\partial Q}{\partial K_{\mathrm{T}}''}\Delta K_{\mathrm{T}}'', \tag{19} \]

where \(\partial P/\partial K_{\mathrm{T}}'\), \(\partial P/\partial K_{\mathrm{T}}''\), \(\partial Q/\partial K_{\mathrm{T}}'\), \(\partial Q/\partial K_{\mathrm{T}}''\) are rectangular matrices, in their structure completely analogous to the submatrices \(\partial P/\partial K_{\mathrm{T}}\), \(\partial Q/\partial K_{\mathrm{T}}\), which respectively are equal to

\[ \frac{\partial P}{\partial K_{\mathrm{T}}'} = \frac{\partial P}{\partial K_{\mathrm{T}}}K_{\mathrm{T}}^{0'}, \quad \frac{\partial P}{\partial K_{\mathrm{T}}''} = \frac{\partial P}{\partial K_{\mathrm{T}}}K_{\mathrm{T}}^{0''}, \quad \frac{\partial Q}{\partial K_{\mathrm{T}}'} = \frac{\partial Q}{\partial K_{\mathrm{T}}}K_{\mathrm{T}}^{0'}, \quad \frac{\partial Q}{\partial K_{\mathrm{T}}''} = \frac{\partial Q}{\partial K_{\mathrm{T}}}K_{\mathrm{T}}^{0''}. \tag{20} \]

The square submatrices \(W_{pU}\), \(W_{qU}\) differ from the matrices \(\partial P/\partial U\), \(\partial Q/\partial U\) by diagonal elements, which are respectively equal to

\[ W_{pUii} = \frac{\partial P_i}{\partial U_i} - \alpha_i, \qquad W_{qUii} = \frac{\partial Q_i}{\partial U_i} - \beta_i. \tag{21} \]

The phase of the voltage vector of one of the nodes \((n)\) of the system must be chosen as the reference, setting for this \(\delta_n = \mathrm{const}\). One of the nodes of the system \((n)\) must be taken as the active-power balance node, whose power increment \(\Delta P_{n\mathrm{s}} = \Delta P_{\mathrm{balance}}\) is forced in nature and, together with \(\Delta U\) and \(\Delta \delta\), is to be determined. Replacing the zero column \(n\) (for the reference voltage vector) in the matrices \(\partial P/\partial \delta\) and \(\partial Q/\partial \delta\) by coefficients with opposite signs for \(\Delta P_{n\mathrm{s}}\) of the disturbance matrices, we obtain the matrices \(W_{p\delta}\) and \(W_{q\delta}\). Correspondingly, in the matrix \(\Delta \delta_{(n)}\), element \(n\) is \(\Delta P_n\). Since

* Equations (17) also make it possible to calculate any initial operating condition for specified node powers, analogously to the calculation of an operating-condition change.

because equation (17) has been compiled taking into account the static characteristics of the elements, the reactive-power operating condition is balanced automatically.

The system of linearized equations (17) can be solved by any known method.

When regulating the voltage only by phase, it is necessary to take into account the additional relation

\[ \Delta K_{Ti}^{\prime\prime} = - \frac{K_{Ti}^{\prime\prime}}{K_{Ti}^{\prime}} \Delta K_{Ti}^{\prime} \tag{22} \]

and

\[ \Delta K_{Ti}^{\prime\prime} = \frac{K_{Ti}^{\prime\prime}}{K_{Ti}^{\prime}} \Delta K_{Ti}^{\prime} \tag{23} \]

— when the transformation ratio is changed only in magnitude.

For transformers with longitudinal voltage regulation, taking into account that \(K_{Ti}^{\prime}/K_{Ti} \simeq 1\) and \(\Delta K_{Ti}^{\prime} \simeq \Delta K_{Ti}\), we obtain the generally accepted simplified expressions:

\[ \frac{\partial P}{\partial K_{Ti}^{\prime}} \Delta K_{Ti}^{\prime} \simeq \frac{\partial P}{\partial K_{Ti}} \Delta K_{Ti}, \qquad \frac{\partial Q}{\partial K_{Ti}^{\prime}} \Delta K_{Ti}^{\prime} \simeq \frac{\partial Q}{\partial K_{Ti}} \Delta K_{Ti}. \tag{24} \]

In those cases where a change (disconnection) of the entire active power of the load concentrated in two-terminal element \(i\) is specified, i.e. \(\Delta P_{i\delta} = -P_i\), in the submatrix \(W_{pU}\) one should set \(a_i = 0\). Similarly, when the entire reactive power of two-terminal element \(i\) is disconnected, the coefficient \(\beta_i\) in the submatrix \(W_{qU}\) is set equal to zero.

In some cases, the change in the voltage magnitude at a node, \(\Delta U_i = \Delta U_{i\delta}\), may be specified as an independent variable. In particular, under astatical voltage regulation \(\Delta U_{i\delta} = 0\), the unknown quantity is the forced component of reactive power \(\Delta Q_{ib}\), which ensures constancy of the voltage. For this it is necessary, in the original equation (17), to make the substitution

\[ \Delta U_i = \frac{1}{\beta_i}\Delta Q_{ib} = \gamma_i \Delta Q_{ib}. \tag{25} \]

In the same way, one can determine the required magnitude of the change in the transformation ratio from a specified value of \(\Delta U_i\) or \(\Delta Q_i\).

When the connection between nodes* \(i\) and \(j\) is opened, it is necessary, using (1) and (2), to calculate the value of the power flowing through the specified connection in the normal operating condition with respect to nodes \(i\) and \(j\). The obtained values, with opposite signs, should be placed in the corresponding positions of the disturbance submatrices; in this case, in the derivatives

\[ \frac{\partial P_{i,j}}{\partial U_{i,j}}, \quad \frac{\partial Q_{i,j}}{\partial U_{i,j}}, \quad \frac{\partial P_{i,j}}{\partial \delta_{i,j}}, \quad \frac{\partial Q_{i,j}}{\partial \delta_{i,j}}, \quad \frac{\partial P_{i,j}}{\partial K_{Tij}}, \quad frac{\partial Q_{i,j}}{\partial K_{Tij}}, \]

\(y_{ij}\) should be set equal to zero. Since in this case the source of the disturbance is not the two-terminal elements \(i\) and \(j\), but the multi-terminal element itself, then \(\eta_i \Delta P_i = 0\), \(\eta_j \Delta P_j = 0\).

Equation (17) remains valid if, in the system under consideration, the generators have ART with respect to external parameters. In particular, if the ART of generator \(i\) is carried out with respect to the mutual angle \(\delta_{in}\), then the unknown will be the increment of active power \(\Delta P_{ib}\), which can be found analogously to the power of the balancing node.

* It is assumed here that the system does not split into independent parts.

When a constraint arises for some parameter $\Pi$, the total increment of the latter must be set equal to zero. The additional equation thus obtained,

\[ \Delta \Pi_i = \Delta \Pi_{iз} + \Delta \Pi_{ib} = 0 \]

together with (17), makes it possible to find the required magnitude of the action $(\Delta \Pi_{iз})$, by means of which the specified operating constraint can be maintained.

After calculating the quantities $\Delta U$ and $\Delta \delta$, one can determine the change in all the other operating parameters and their new values. For small changes in normal operation, the slope coefficients $\alpha, \beta, \eta$ change only slightly and may be taken as constant in the calculation. Otherwise, the values of the coefficients must be corrected as the operating condition changes.

Equivalencing parts of an electric system by operating parameters

The order of equation (17) can be reduced considerably if individual parts of the electric system are represented by an operating equivalent—generalized static response characteristics [5, 9]. Such equivalencing is especially convenient when describing the system operating condition by equations (17) and when calculating limiting operating conditions, according to static-stability criteria, for chain networks (for example,

Fig. 2. Equivalencing a system section by a two-terminal network

long-distance transmission lines with intermediate systems), since it makes it possible, in generalized form, to take into account the response of individual parts of the system and to isolate directly the object of disturbance or control.

For the equivalent parts of the system, represented by a two-terminal network, it is expedient to use the same form of notation for the static characteristics with respect to active power as for a load (10), and with respect to reactive power as for a generator (12), i.e.,

\[ Q_{ic} = Q(U_i, P_{ic}), \qquad P_{ic} = P(U_i). \tag{26} \]

The possibility of such a replacement of an electric power system (or of its part) follows from the fact that the equivalent part, in the general case, represen-

is a combination of generator and load buses and of the passive part—the electrical network—which is the transmission link between the indicated elements and the equivalencing bus.

The slope coefficients of the generalized static characteristics of the part of the system being equivalenced (26) are determined for a fixed network diagram, AVR settings, specified slope coefficients of the static characteristics of generators, synchronous compensators, and loads, and a known initial operating condition. In addition, the balancing bus must be specified.

Let us consider the calculation of the quantities \(\alpha_{ic}\), \(\beta_{ic}\), \(\eta_{ic}\).

The coefficients \(\alpha_{ic}\) and \(\beta_{ic}\) characterize the response of the system, respectively in active and reactive power, to a change in voltage at bus \(i\) by an external source (Fig. 2).

The matrix of coefficients \(\alpha_{ic}\) for all nodal points of the power system can be obtained from equation (17), assuming that a disturbance \(\Delta P_{ic}\) is applied successively to the buses; in this case the phase of the voltage vector must be fixed at the equivalencing bus, which is at the same time the bus balancing the reactive power, i.e.,

\[ \alpha_{ic}=\frac{dP_{ic}}{dU_i}=\frac{D_\alpha}{A_{i,i}}, \tag{27} \]

where \(D_\alpha\) is the determinant, and \(A\) is the cofactor of the matrix \(W_\alpha\)

\[ W=\left\|\begin{matrix} W_{pu} & W_{p\delta}\\ W_{qu} & W_{q\delta} \end{matrix}\right\|. \tag{28} \]

Since the source of the disturbance is the “imbalance” applied to bus \(i\) with the two-terminal element \(i\) disconnected, the matrix \(W_\alpha\) differs from \(W\) in that in it \(\alpha_i\), \(\beta_i\) are equal to 0. In the matrix \(W_\alpha\), all elements of column \(n+i\) (the coefficients at the balancing power \(\Delta Q_{ic}\)) are equal to zero, except for \(a_{n+i,n+i}=-1\).

Similarly, assuming that a disturbance \(\Delta Q_{ic}\) is applied successively to the buses of the system, one can calculate the matrix of coefficients

\[ \beta_{ic}=\frac{\partial Q_{ic}}{\partial U_i}=\frac{D_\beta}{A_{i,n+i}}. \tag{29} \]

The coefficient \(\eta_{ic}\) shows by how many units it is necessary to change the reactive power transmitted to bus \(i\) by the system so that, when its active power changes by one unit, the voltage \(U_i\) remains constant. The value of the coefficient \(\eta_{ic}\) depends on the distribution of the increment \(\Delta P_i\) among the generators and loads of the power system being equivalenced. Therefore, when determining \(\eta_{ic}\), the law (proportions) of distribution of the power increment \(\Delta P_i\) among the buses of the power system must be specified. This increment may be distributed among several generation and load buses, or may be assigned entirely to the balance bus \((n)\). In the latter case, the value \(\eta_{ic}\) can be found from (17), assuming \(\Delta U_i=0\) and \(\Delta \delta_i=0\):

\[ \eta_{ic}=\frac{\partial Q_{ic}}{\partial P_{ic}}\bigg|_{(U_i=\mathrm{const})} =\frac{A_{i,i}}{D_\eta} =-\frac{A_{i,i}}{A_{n+i,n}+\eta_n A_{n+i,2n}}. \tag{30} \]

Let us note that, in general, in the slope coefficient \(\beta_{ic}\) (29) it is possible, in “pure” form, to single out the components caused by the influence of the regulating properties of the active and reactive powers of loads (or of a bus separately) and to estimate their magnitude.

Calculation of the Static Stability of the System

When the operating state of an electric power system changes, it becomes necessary to check the stability of the new steady state, and also to determine the parameters of the limiting state and the corresponding stability margins [8].

The stability of a stationary state of an electric power system without allowance for self-swinging can be judged by the criterion of positivity of the constant term of the characteristic equation [9], or by practical criteria [6].

The initial matrix equation (17) makes it possible, from the known parameters of the steady state, to calculate the value of any practical criterion.

Assuming that the source of the reactive-power disturbance is two-terminal element \(i\), analogously to (29) we have

\[ \frac{dQ_i}{dU_j}=\frac{D}{A_{j,n+i}} . \tag{31} \]

Fixing the phase of the voltage vector of the balancing node \((n)\), in an analogous manner we obtain the derivative of the power of node \(P_i\) with respect to the mutual angle \(\delta_{in}\)

\[ \frac{dP_i}{d\delta_{in}}= \frac{D}{A_{n+i,i}+\eta_i A_{n+i,n+i}} . \tag{32} \]

In expressions (31), (32), \(D\) is the determinant, and \(A\) is the algebraic complement of the matrix \(W\) (28).

The limiting state of the system corresponds to equality to zero of the practical criteria [6], or, equivalently, \(D=0\).

The determinant \(D\) of the matrix \(W\) changes its sign only upon crossing the boundary from a stable state of the system into an unstable one, or conversely. Therefore, if the sign of the determinant \(D\) of equation (17) is known for a definitely stable state, then there is no need to check the new steady state by expressions (31) or (32), since the presence or absence of static stability can be judged directly in the process of finding the state from the sign of the determinant \(D\).

For a complex electric system, under the conditions of aperiodic stability, there corresponds a set of limiting states, which are determined by: 1) the corresponding values of the parameters of the adopted initial state; 2) the choice of the station to be checked for stability, or of the section of the system for determining limiting power flows; 3) the choice of the balancing node (station) in the system.

The limiting state of the system is found by successively increasing the loading of the normal initial state by redistributing the active powers of the generating stations; in this case the active powers of \((S-2)\) stations remain unchanged, one is loaded, and another is the balancing station [9]. The state can be made heavier in accordance with the conditions set forth above by changing (disturbing) various system parameters.

1. By a discrete increment of the active power of the loaded station \(\Delta P_{f3}\).

When solving equation (17), the mismatch of node \(i\) in active and reactive power at the \(k\)-th iteration within each loading step is determined:

\[ \text{for } i=1,2,\ldots,n-1\quad (i\ne f) \]

\[ \left. \begin{aligned} F_{pi} &= P_i + \sum_{\kappa=1}^{\kappa-1} \alpha_i \Delta U_i - U_i^2 g_{ii} - U_i \sum_{\substack{j=1\\ j\ne i}}^{n} U_j y_{ij}\sin(\hat{\delta}_i-\hat{\delta}_j-\alpha_{ij}),\\ F_{qi} &= Q_i + \sum_{\kappa=1}^{\kappa-1} \beta_i \Delta U_i - U_i^2 b_{ii} + U_i \sum_{\substack{j=1\\ j\ne i}}^{n} U_j y_{ij}\cos(\hat{\delta}_i-\hat{\delta}_j-\alpha_{ij}), \end{aligned} \right\} \tag{33} \]

for \(i=f, n\) (\(f\) is the number of the loaded station, \(n\) the balancing one)

\[ \left. \begin{aligned} F_{pi} &= P_i + \Delta P_i - U_i^2 g_{ii} - U_i \sum_{\substack{j=1\\ j\ne i}}^{n} U_j y_{ij}\sin(\hat{\delta}_i-\hat{\delta}_j-\alpha_{ij}),\\ F_{qi} &= Q_i + \gamma_i F_{pi} + \sum_{\kappa=1}^{\kappa-1} \beta_i \Delta U_i - U_i^2 b_{ii} + U_i \sum_{\substack{j=1\\ j\ne i}}^{n} U_j y_{ij}\cos(\hat{\delta}_i-\hat{\delta}_j-\alpha_{ij}). \end{aligned} \right\} \tag{34} \]

As noted above, according to practical criteria, the parameters of the limiting operating condition of the system under the conditions of aperiodic stability correspond to \(\det W=0\). At the same time, this same condition is the criterion for instability of convergence of the Newton iteration process [10]*. Consequently, on the one hand, with the adopted method of loading the operating condition, the Newton method cannot be used to calculate the parameters of the limiting operating condition of the system under study (although an operating condition sufficiently close to the limiting one can be calculated); on the other hand, violation of convergence of the iterative calculation of equation (17) corresponds in this case to the limit of static stability of the system**.

The indicated difficulties associated with finding the limiting operating conditions of the system can be eliminated without changing the method of solving the equations adopted as the basis in the preceding discussion. For this it is necessary to pass to another way of loading the operating condition (to change the disturbing parameter).

1. Let us take as the independent variable a prescribed increase in the angle of displacement between the loaded and balancing stations, \(\delta_{fn}\) (\(\delta_n=\mathrm{const}\)); in this case the increment of active power of node \(f\) is of a forced character (\(\Delta P_{fb}\)). Accordingly, the matrix \(W\) in (17) must be replaced by the matrix \(W_{\delta}\), which differs from the former in that the elements of column \(n+f\) are the coefficients at \(\Delta P_f\) (all elements of the indicated column are equal to zero except \(a_{n+f,f}=-1,\ a_{n+f,n+f}=-\beta_f\)).

Since \(\det W \ne \det W_{\delta}\), for the parameters of the limiting operating condition of the system the Newton method provides quite stable solutions of the equations for the entire manifold of steady-state operating conditions of the system, which are both statically stable and unstable.

Within each step of loading the operating condition, the first iteration is calculated with the vector function of disturbance equal to

\[ \begin{Vmatrix} G_p\\ G_q \end{Vmatrix} = - \begin{Vmatrix} \dfrac{\partial P}{\partial \delta_f}\\[4pt] \dfrac{\partial Q}{\partial \delta_f} \end{Vmatrix} \Delta \delta_{f3}, \tag{35} \]

* Strictly speaking, \(\det W \ne 0\), since the parameters of the limiting operating condition have not been found; however, the closeness of the matrix \(W\) to a singular one is the reason why another condition for convergence of the process is violated: at some stage of the calculation the absolute value of the difference between the values of the sought root and the approximate root of the \(k+1\)-st approximation is greater than the corresponding value of the \(k\)-th approximation, i.e., the process diverges.

** We note that the method of steepest descent in this case has stable convergence.

for subsequent iterations, the mismatch for each of the buses is determined from the corresponding expressions (33), (34).

We note that the operating condition can be made heavier by a specified change (decrease) of the voltage at one of the buses of the system, analogously to the case considered above.

Conclusions

1. The equations obtained make it possible to find the magnitude of the control action that ensures the desired change in the operating condition, and to evaluate the response of the system, by means of transfer functions, to a disturbance of any operating parameter.

2. The static stability of the system in the absence of self-excitation can be judged in the process of calculating the operating condition.

3. In the general case, a violation of convergence of the iterative calculation process is not a criterion of static instability of the system. The interrelation between static instability (stability) of the system and disruption of the convergence process in solving the equations should be established only through the interrelation of their criteria, which depends on the way in which information in the electric system is specified, on the method used for solving the equations, and on the choice of independent parameters (arguments) of the operating condition. Therefore, by varying the “way,” “method,” and “choice,” it is always possible to obtain a stable solution for the entire variety of system operating conditions.

Appendix

Formulas for partial derivatives

1. Derivatives of the power of bus \(i\) with respect to the magnitude and phase of the voltage of this bus

\[ \frac{\partial P_i}{\partial U_i}=2U_i g_{ii}+A_i,\qquad \frac{\partial Q_i}{\partial U_i}=2U_i b_{ii}-B_i, \tag{п-1} \]

\[ \frac{\partial P_i}{\partial \delta_i}=U_i B_i,\qquad \frac{\partial Q_i}{\partial \delta_i}=U_i A_i, \tag{п-2} \]

where

\[ g_{ii}=\sum_{\substack{j=1\\ j\ne i}}^{n} y_{ij}\sin \alpha_{ij}K_{Tij}^{2},\qquad b_{ii}=\sum_{\substack{j=1\\ j\ne i}}^{n} y_{ij}\cos \alpha_{ij}K_{Tij}^{2}, \tag{п-3} \]

\[ A_i=\sum_{\substack{j=1\\ j\ne i}}^{n} U_j y_{ij}\sin(\delta_i-\delta_j-\alpha_{ij})K_{Tij}, \]

\[ B_i=\sum_{\substack{j=1\\ j\ne i}}^{n} U_j y_{ij}\cos(\delta_i-\delta_j-\alpha_{ij})K_{Tij}. \tag{п-4} \]

If the branch admittance is connected to the bus directly (and not through a transformer), then the corresponding value of \(K_T\) in (п-3) is equal to unity.

1. The derivatives of the power of node \(i\) with respect to the magnitude and phase of the voltage of adjacent nodes \(j = 1, 2, \ldots, n\) \((j \ne i)\)

\[ \frac{\partial P_i}{\partial U_j}=c_i,\quad \frac{\partial Q_i}{\partial U_j}=D_i,\quad \frac{\partial P_i}{\partial \delta_j}=U_jD_i,\quad \frac{\partial Q_i}{\partial \delta_j}=-U_jC_i, \tag{п-5} \]

where

\[ C_i=U_i y_{ij}\sin(\delta_i-\delta_j-\alpha_{ij})K_{Tij},\quad D_i=-U_i y_{ij}\cos(\delta_i-\delta_j-\alpha_{ij})K_{Tij}. \tag{п-6} \]

1. The derivatives of power with respect to the transformation ratio \(K_{Tif}\) of branch \(if\); the branch admittance is referred to node \(i\) in the equivalent circuit, and the transformer is directly connected to node \(f\) (Fig. 1):

a) for node \(i\)

\[ \frac{\partial P_i}{\partial K_{Tif}} = U_iU_fy_{if}\sin(\delta_i-\delta_f-\alpha_{if}), \]

\[ \frac{\partial Q_i}{\partial K_{Tif}} = -\,U_iU_fy_{if}\cos(\delta_i-\delta_f-\alpha_{if}); \tag{п-7} \]

b) for node \(f\)

\[ \frac{\partial P_f}{\partial K_{Tif}} = 2U_f^2g_{if}K_{Tif} + U_fU_i y_{if}\sin(\delta_f-\delta_i-\alpha_{fi}), \tag{п-8} \]

\[ \frac{\partial Q_f}{\partial K_{Tif}} = 2U_f^2b_{if}K_{Tif} - U_fU_i y_{if}\cos(\delta_f-\delta_i-\alpha_{fi}). \tag{п-9} \]

The transformation ratios are everywhere determined in the direction from the node to which the branch admittance is connected to the node to which the transformer is connected (Fig. 1)

\[ K_{Tif}=U_i/U_f. \]

Correspondingly, the referred values of the voltages and admittances of the nodes in (1), (2) are equal to

\[ U_{f(i)}=U_fK_{Tif},\quad U_{i(f)}=U_i/K_{Tif},\quad y_{fi}=y_{if}K_{Tif}^{2},\ \text{etc.} \]

REFERENCES

1. N. A. Melnikov, L. A. Soldatkina. Matrix Methods for Calculating Operating Conditions of Closed Electric Networks. GEI, Moscow, 1965.

2. V. A. Venikov, N. D. Anisimova, N. K. Krug. On the question of determining self and mutual admittances when varying the parameters of a part of an equivalent circuit or combining two systems into one. Elektrichestvo, No. 1, 1970.

3. V. A. Venikov, N. D. Anisimova, R. Mitske, T. Shchechina. Investigation by means of a digital computer of the static stability of complex automatically regulated electric systems. Elektrichestvo, No. 11, 1967.

4. N. A. Melnikov. Methods for calculating operating conditions for circuits containing transformation elements with complex parameters. Izv. AN SSSR. Energetika i transport, No. 4, 1964.

5. L. M. Zeisman. Some questions of controlling the steady-state conditions of thermal power plants with respect to voltage and reactive power. Author’s abstract of candidate’s dissertation, Moscow, 1967.

6. I. M. Markovich. Operating Conditions of Power Systems. GEI, Moscow, 1963.

7. N. A. Melnikov. A Matrix Method for the Analysis of Electric Networks. “Energiya,” Moscow, 1966.

8. D. I. Azar’ev, V. A. Venikov et al. Basic principles for determining the stability of power systems. Elektrichestvo, No. 11, 1963.

9. P. S. Zhdanov. Stability of Electric Systems. Gosenergoizdat, Moscow, 1948.

10. B. P. Demidovich, I. M. Maron. Fundamentals of Computational Mathematics. Fizmatgiz, Moscow, 1963.