COMPUTERIZED DEVICE FOR THE MARGINAL STATE ASSESSMENT OF POWER SYSTEMS

Abstract

The required technical result aimed to raise a speed of device response is achieved in the device, which includes a group of operative memory blocks, a data acquisition block, memory block and marginal state assessment block, where operative memory blocks group outputs are connected with data acquisition block inputs, which output is connected, with marginal state assessment block input, which output is connected with the memory block input, and first and second outputs of marginal state assessment block are connected, respectively, with first and second inputs of the memory block.

Claims

1. A computerized device for a marginal state assessment of power systems, comprising: a group of operative memory blocks, a data acquisition block, an output memory block, and a marginal state assessment block, wherein an output of the operative memory blocks of the group are connected with an input of the data acquisition block, wherein an output of the data acquisition block is connected with an input of the marginal state assessment block, wherein a first output of the marginal state assessment block is connected with an input of the output memory block, wherein the marginal state assessment block comprises a series of calculators, said series of calculators comprising a calculator of a Lagrange multiplier vector, a calculator of loading factor ultimate increment, a calculator of sign of power flow matrix determinant, a calculator of increments and corrections, and a calculator of convergence check, wherein an input of the calculator of a Lagrange multiplier vector is the input of the marginal state assessment block, wherein an output of the calculator of a Lagrange multiplier vector is an input of the calculator of loading factor ultimate increment, an output of the calculator of loading factor ultimate increment is an input of the calculator of sign of power flow matrix determinant, an output of the calculator of sign of power flow matrix determinant is an input of the calculator of increments and corrections, and an output of the calculator of increments and corrections is an input of the calculator of convergence check, wherein a first output of the calculator of convergence check is the first output of the marginal state assessment block and a second output of the calculator of convergence check is connected with an additional input of the calculator of a loading factor ultimate increment and further connected to the calculator of a sign of power flow matrix determinant, further connected to a calculator of bifurcation, wherein an output of the calculator of bifurcation is a second output of the marginal state assessment block, wherein the first and the second outputs of the marginal state assessment block are connected correspondingly with the first and the second inputs of the output memory block.

2. A method for performing a marginal state assessment of a power system, comprising: receiving one or more outputs from a group of operative memory blocks, sending the one or more outputs to an input of a data acquisition block, sending an output of the data acquisition block to an input of a marginal state assessment block, within said marginal assessment block: first calculating a Lagrange multiplier vector using the output of the data acquisition block, second calculating a loading factor ultimate increment using the output of the Lagrange multiplier vector, third calculating a sign of power flow matrix determinant using the loading factor ultimate increment, fourth calculating increments and corrections, fifth calculating a convergence check, thus forming a first output of the marginal state assessment block, said first output of the marginal state assessment block being a first input of an output memory block, additionally sending said first output of the marginal state assessment block to calculate a second loading factor ultimate increment, said loading factor ultimate increment being used to calculate a second sign of power flow matrix determinant, said second sign of power flow matrix determinant being used to calculate a bifurcation, wherein the bifurcation is a second output of the marginal state assessment block, said second output of the marginal state assessment block being a second input of the output memory block, and making a final determination within the output memory block, the final determination comprising an indicator of at least one marginal state of the power system.

3. The device of claim 1, wherein a constraint violation determined by the marginal state assessment block causes a switch in a corresponding generator.

4. The device of claim 1, wherein an output of the marginal state assessment block causes a switch in a corresponding generator bus from PV to PQ.sup.+.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0015] The following is presented on the figures:

[0016] On FIG. 1—functional scheme of computerized device for the marginal state assessment of power systems;

[0017] On FIG. 2—table of results of marginal states assessment for an example of two-node system.

DETAILED DESCRIPTION OF THE INVENTION

[0018] The computerized device for the marginal state assessment of power systems (FIG. 1) includes a group of 1-1 . . . 1-n operative memory blocks, a data acquisition block 2, a marginal state assessment block 3 and an output memory block 4, in addition the outputs of 1-1 . . . 1-n operative memory blocks of the group are connected with the input of data acquisition block 2, the output of this block 2 is connected with the input of a marginal state assessment block 3, the first and the second outputs of this block 3 connected correspondingly with the first and the second inputs of an output memory block 4.

[0019] In the computerized device for the marginal state assessment of power systems a marginal state assessment block is made in a form of connected in series a calculator 5 of Lagrange multiplier vector with the input of this calculator being the input of a marginal state assessment block 3, a calculator 6 of loading factor ultimate increment, a calculator 7 of power flow matrix determinant sign, a calculator 8 of increments and corrections and a calculator 9 of convergence check, the first output of this calculator is the first output of a marginal state assessment block 3 and the second output of this calculator is connected with an additional input of a calculator 6 of loading factor ultimate increment, and also a calculator 10 of bifurcation the input of this calculator is connected with a additional output of a calculator 7 of power flow matrix determinant sign and the output of it is the second output of a marginal state assessment block 3, besides the first and the second outputs of a marginal state assessment block 3 are connected correspondingly with the first and the second inputs of an output memory block 4.

[0020] The proposed technical decision includes the elements characterized on the functional level and a described form of realization is supposed to use a programmable (adjustable) multifunctional facility, so further on description of its operation the information is presented that confirms an opportunity of this facility to perform the specific function required by the given technical decision particularly to solve the algorithms and mathematic expressions.

[0021] The computerized device for the marginal state assessment of power systems operates in the following way.

[0022] Preliminary let us perform a theoretical justification of the device operation algorithm.

[0023] The assessment of marginal states (MS) plays a key role in the analysis, planning and control of power systems. The MS is such a stable steady state for which an arbitrarily small change in any of operating quantities in an unfavorable direction causes a voltage collapse or a loss of synchronism by units.

[0024] One of the most commonly used approaches of the MS search consists of sequential load flow solutions with the given step along the specified trajectory until a divergence of load flow solution is obtained, which in turn is corrected by the binary search or by the damped Newton method. However, as soon as the network comes close to the condition of instability, a divergence in load flow solution may occur that is caused by the numerical problems of ill-conditioned Jacobian. That is why, in continuation methods the load flow parameterization, or the normalized iteration step change, are widespread. In this case the Jacobian of intermediate states and MS is not singular.

[0025] The MS model of power system can be presented as follows

min−t  (1)

subject to

ΔF(X,Y+tdY)=0.  (1a)

[0026] Here t is a scalar variable typically referred to as the loading factor.

[0027] System (1a) represents the active power balance equations at the PV and the PQ buses, and the reactive power balance equations at the PQ buses; the vector X represents the bus voltage angles (except the slack bus angle), and the PQ bus voltage magnitudes; the vector Y represents the active and reactive power at each load bus and the active power generated at each generator bus; the direction vector dY represents the user-defined changes in active and reactive power demand and the changes in active power generation.

[0028] The Lagrange function for (1)-(1a) can be presented as

L=−t+ΔF(X,Y+tdY).sup.Tλ,   a. (2)

[0029] where λ is a vector of auxiliary solution variables called Lagrange multipliers.

[0030] Differentiation of this function with respect to all variables yields

∇.sub.XL=[J].sup.Tλ=0;  (3)

∇.sub.λL=ΔF(X,Y+tdY)=0;  (3a)

∇.sub.tL=−1+dY.sup.Tλ=0.  (3b)

[0031] Here ∇.sub.XL=[∂L/∂X].sup.T is a gradient of the Lagrange function with respect to vector X, and [J]=[∂ΔF/∂X] is the load flow Jacobian.

[0032] Condition (3b) ensures that the vector λ will not be equal to zero at a solution point of (1)-(1a). That is why condition (3) determines the load flow Jacobian singularity, i.e., an MS.

[0033] In the MS the load flow Jacobian is singular. Contrary to it the matrix of linearized equations (3)-(3b)

[00001]$\begin{array}{cc}\left[\begin{array}{ccc}H& J^T& 0\\ J& 0& \mathrm{dY}\\ 0& {\mathrm{dY}}^T& 0\end{array}\right]\left[\begin{array}{c}\Delta .Math..Math.X\\ \mathrm{\Delta \lambda }\\ \Delta .Math..Math.t\end{array}\right]=-\left[\begin{array}{c}{\nabla }_X.Math.L\\ {\nabla }_{\lambda }.Math.L\\ {\nabla }_t.Math.L\end{array}\right]& \left(4\right)\end{array}$

[0034] is not singular at the solution point of (1)-(1a). Here [H]=└∂.sup.2ΔF.sup.Tλ/∂X.sup.2┘ is a matrix of second partial derivatives (Hessian).

[0035] Equations (3)-(3a) can be obtained directly from the implicit function theorem. In this case they have got the name of the point of the collapse (PoC) method, and in Russia they are called the MS equations. In existing PoC modifications the condition (3) in some cases is written for the right eigenvector of the Jacobian, which corresponds to the zero eigenvalue, while the Euclidean or the max norm of vector λ is used as the addend in (3b).

[0036] Researches have shown that the good initial guesses for the system variables, particularly the eigenvectors, are essential, otherwise the Newton approach for obtaining the solution to the PoC equations (3-3b) either yields undesirable results or does not converge. That is why, in spite of the fact that the PoC method surpasses the continuation methods in case of the successful solution, because of unreliability of the PoC method, many researchers prefer to use continuation methods, at least, for finding an initial MS.

[0037] Consider a 2-bus system with PV-bus k and slack bus m. In this case the problem (1)-(1a) will look as follows:

min−t

subject to

ΔP.sub.k=P.sub.k+tdP.sub.k−V.sub.k.sup.2|Y.sub.km|sin α.sub.km−V.sub.kV.sub.m|Y.sub.km|sin(δ.sub.km−α.sub.km)=0,

[0038] where t is the loading factor; δ.sub.k, V.sub.k are the voltage angle and magnitude at bus k; P.sub.k is the active power at bus k; dP.sub.k is the given “direction” of active power change at bus k; |Y.sub.km| is the series admittance magnitude between buses k and m; α.sub.km=∠Y.sub.km+π/2 is a loss angle; δ.sub.km=δ.sub.k−δ.sub.m.

[0039] The Lagrange function for this problem can be presented as

L=−t+[P.sub.k+tdP.sub.k−V.sub.k.sup.2|Y.sub.km|sin α.sub.km−V.sub.kV.sub.m|Y.sub.km|sin(δ.sub.km−α.sub.km)]λ.

[0040] The first-order optimality conditions are

∇.sub.δ.sub.kL=−V.sub.kV.sub.m|Y.sub.km|cos(δ.sub.km−α.sub.km)λ=0;  (5)

∇.sub.λL=ΔP.sub.k=P.sub.k+tdP.sub.k−V.sub.k.sup.2|Y.sub.km|sin α.sub.km−V.sub.kV.sub.m|Y.sub.km|sin(δ.sub.km−α.sub.km)=0;  (5a)

∇.sub.tL=−1+dP.sub.kλ=0.  (5b)

[0041] Equation (5b) directly determines the Lagrange multiplier λ=1/dP.sub.k. Therefore if to use it as the initial guess of multiplier λ then only two linearized equations will be really used during iterations.

[V.sub.kV.sub.m|Y.sub.km|sin(δ.sub.km−α.sub.km)λ]Δ.sub.k=V.sub.kV.sub.m|Y.sub.km|cos(δ.sub.km−α.sub.km)λ;  (6)

[−V.sub.kV.sub.m|Y.sub.km|cos(δ.sub.km−α.sub.km)]Δδ.sub.k+dP.sub.kΔt=−ΔP.sub.k.  (6a)

[0042] According to (6)

Δδ.sub.k=ctg(δ.sub.km−α.sub.km)  (7)

[0043] Substitution of (7) in (6a) yields

Δt=−ctg(δ.sub.km−α.sub.km)∇.sub.δ.sub.kL−ΔP.sub.k/dP.sub.k.  (8)

[0044] Depending on the initial state and the given direction dP.sub.k the multiplier in Δδ.sub.k (6), which corresponds to the Hessian [H] in (4), can be positive or negative. The change of the state (7) will go only in one direction irrespective of the given direction dP.sub.k. In the case of δ.sub.km>α.sub.km the voltage angle δ.sub.km will increase, hence the MS will correspond to the maximum generation of active power at bus k. In the case of δ.sub.km<α.sub.km the voltage angle δ.sub.km will decrease, and the MS will correspond to the maximum load at bus k. The steady state with δ.sub.km=α.sub.km is the “border-state” (BS) between the attraction regions to these MS. The determinant of matrix of equations (6)-(6a) is equal to V.sub.kV.sub.m|Y.sub.km|sin(δ.sub.km−α.sub.km). Therefore in the BS i.e., when δ.sub.km=α.sub.km, the matrix of linearized equations (6)-(6a) is singular. In the BS neighbourhood this matrix is ill-conditioned. Therefore in the case of the “light” initial state, where δ.sub.km≈α.sub.km, an increase in Δδ.sub.k (7) and Δt (8) can be abnormal. As opposed to it in the MS the determinant has the highest absolute value. In the MS proximity the matrix of system (6)-(6a) is well-conditioned and the Newton method has a quadratic convergence.

[0045] In rectangular coordinates the MS will also depend on the initial steady state. An angle shift for more than one period is impossible, but the loading factor t and the bus voltage can essentially exceed required values. In the MS proximity the convergence speed turned out lower than in polar coordinates. It can be supposed that the usage of t.sup.2 dP is a more successful alternative than tdP. However, it does not yield the desired result. If the initial Hessian [H] is not positive definite, the oscillatory iterative process emerges around the BS δ.sub.km=α.sub.km. In the case of the positive definite matrix the convergence speed turned out noticeably worse than with tdP.

[0046] Thus, the presented analysis shows that for realization of the problem (1)-(1a) it is necessary to solve two problems. Firstly, the iterative process should generate a change of variables in the “correct” direction. Secondly, using the system of linear equations (4), it is necessary to apply mechanisms which take into account ill-conditioned matrix of this system. The theory of nonlinear programming has many procedures for these purposes. However, a common realization of these procedures has not yet allowed to create a computationally efficient algorithm of solution of the constrained minimization problem (1)-(1a).

[0047] There are two classical approaches for globalizing the locally convergent algorithm: line search procedures, or use of trust regions. The first approach requires the usage of a merit function. The selection of the merit function is ambiguous and strongly influences the speed of algorithm convergence. The second approach requires computationally costly procedures to support, e.g., positive definiteness of the Hessian, or in order that inertia of the matrix (4) satisfied the given conditions.

[0048] At the same time, the use of peculiarities of the system structure of linear equations (4) and the second-order optimality conditions for problem (1)-(1a) allows to create the simple, fast and reliable method of searching MS in a given direction of power change.

[0049] Let us consider a method of MS determination with a single slack bus.

[0050] Present the linearized equations system (4) in the following form

[H]ΔX+[J].sup.TΔλ=−∇.sub.XL;  (9)

[J]ΔX+dYΔt=−ΔF;  (9a)

dY.sup.TΔλ=−∇.sub.tL.  (9b)

[0051] The system of nonlinear equations (3)-(3b) is the first-order optimality condition for problem (1)-(1a). According to the theory of nonlinear programming, a solution of the problems (1-1a) will be a point of the strict local minimum, if

Z.sup.T[H]Z>0  (10)

[0052] for any vector Z which satisfies the equality constraint

[J]Z+dY.Math.Z.sub.t=0,  (10a)

[0053] where Z.sub.t is any scalar. Requirements (3)-(3b) and (10)-(10a) are known as the second-order sufficiency conditions.

[0054] Unlike an unconstrained minimization problem, conditions (10)-(10a) do not require the positive definite Hessian [H] at the solution point. It is enough to have the positive definite Hessian [H] along the direction Z.

[0055] Comparing (10a) with (9a) shows that (10a) is system (9a) with ΔF=0 i.e., when (9a) is homogeneous. Therefore, in order to take into account condition (10) during the problem solution (1)-(1a), represent vector ΔX as the sum of two vectors

ΔX=ΔX.sub.dYΔt+ΔX.sub.ΔF,  (11)

[0056] where ΔX.sub.dY and ΔX.sub.ΔF are solution vectors of the following systems of linear equations:

[J]ΔX.sub.dY=−dY;  (12)

[J]ΔX.sub.ΔF=−ΔF.  (12a)

[0057] Substitution of (11) in (9) yields

Δλ=−[J.sup.T].sup.−1(∇.sub.XL+[H](ΔX.sub.dYΔt+ΔX.sub.ΔF)).  (13)

[0058] In turn, substitution of (13) in (9b) yields

Δt=−(α+ΔX.sub.dY.sup.T[H]ΔX.sub.ΔF)/(ΔX.sub.dY.sup.T[H]ΔX.sub.dY),  (14)

where a=∇.sub.tL+∇.sub.XL.sup.TΔX.sub.dY.

[0059] Comparing (10a) with (12) shows that ΔX.sub.dY=Z/Z.sub.t. Hence ΔX.sub.dY.sup.T[H]ΔX.sub.dY=Z.sup.T[H]Z/Z.sub.t.sup.2. Therefore, to provide a change of variables in the “right” direction, according to requirement (10) it is necessary that the denominator in (14) was positive. If it is not so, it is necessary to “correct” it. The simplest technique that at the same time has a theoretical substantiation is the increase of the diagonal elements of the Hessian [H] by a positive β>0 in (9). In this case (13) is changed and (14) is transformed to

[00002]$\begin{array}{cc}\Delta .Math..Math.t=-\frac{a+\Delta .Math..Math.X_{\mathrm{dY}}^T\left[H\right].Math.\Delta .Math..Math.X_{\Delta .Math..Math.F}+\mathrm{\beta \Delta }.Math..Math.X_{\mathrm{dY}}^T.Math.\Delta .Math..Math.X_{\Delta .Math..Math.F}}{\Delta .Math..Math.X_{\mathrm{dY}}^T\left[H\right].Math.\Delta .Math..Math.X_{\mathrm{dY}}+\beta .Math.{.Math.\Delta .Math..Math.X_{\mathrm{dY}}.Math.}_2^2},& \left(14.Math.a\right)\end{array}$

[0060] where ∥•∥.sub.2 is the Euclidean norm.

[0061] The β value can be selected by the different ways, the basic requirement is

ΔX.sub.dY.sup.T[H]ΔX.sub.dY+β∥ΔX.sub.dY∥.sub.2.sup.2>0.  (15)

[0062] But this requirement is insufficient for reliable convergence of the method.

[0063] The analysis of the 2-bus system has shown that the matrix of linear equations (4), so (9)-(9b), is singular in the BS. For (14) it means that the denominator will be equal to zero. In the BS neighborhood this matrix is ill-conditioned and the denominator (14) will be much less than 1. As a consequence, a change of or (14) and also ΔX (11) and Δλ (13) will be excessive. According to (14a) β allows to reduce Δt, i.e., to improve the condition number of system (9)-(9b). For this purpose it is possible to take advantage of the trust regions ideology. If one way or another a permissible step size Δt.sub.max is specified, then fl can be obtained directly from (14a)

[00003]$\begin{array}{cc}\beta =\max .Math.\left\{0,-\frac{a+\Delta .Math..Math.X_{\mathrm{dY}}^T\left[H\right].Math.\Delta .Math..Math.X_{\Delta .Math..Math.F}+\Delta .Math..Math.X_{\mathrm{dY}}^T\left[H\right].Math.\Delta .Math..Math.X_{\mathrm{dY}}.Math.\Delta .Math..Math.t_{\max }}{\Delta .Math..Math.X_{\mathrm{dY}}^T.Math.\Delta .Math..Math.X_{\Delta .Math..Math.F}+{.Math.\Delta .Math..Math.X_{\mathrm{dY}}.Math.}_2^2.Math.\Delta .Math..Math.t_{\max }}\right\}.& \left(16\right)\end{array}$

[0064] Since the distance to the MS is not known beforehand and the dY size can be chosen arbitrarily, the best solution would be to reassign the Δt.sub.max value adaptively after the calculation of ΔX.sub.dY from (12) in order to limit the maximum change of the voltage angles or other operating parameters at the iteration.

[0065] After determination of Δt from (14a), ΔX is obtained from (11), and Δλ—from the linear equations system solution

[J].sup.TΔλ=−(∇.sub.XL+[H]ΔX+βΔX).  (17)

[0066] The loading factor t is linearly included in (3a). After iteration it will change according to linearized expressions (9a) where nonlinear changes are not considered. For this reason, after each iteration the value of t can be corrected by

Δt.sub.c=ΔF.sup.TdY/∥dY∥.sub.2.sup.2.  (18)

[0067] It reduces the Euclidean norm of the power mismatch vector and this vector becomes orthogonal to dY, that also improves the condition number.

[0068] Let's consider a method of MS determination with a distributed slack bus. In this case the load flow equations (1a) take the following form:

ΔF(X,Y+P.sup.Sα+tdY)=0,  (19)

[0069] where P.sup.S is the slack bus power; α represents the bus participation factors in the distributed slack bus with α.sub.k≧0 and Σα.sub.k=1.

[0070] When (19) is implemented for the problem (1)-(1a), (3b) will remain unchanged, (3a) will be replaced by (19), and (3) will be complemented by the following equation:

∇.sub.P.sub.sL=λ.sup.Tα=0

[0071] Vectors ΔX.sub.dY and ΔX.sub.ΔF used in (11), (14a), and (16) are obtained by solving the following linear equations:

[00004]$\left[J_{\mathrm{LF}}^S\right]\left[\begin{array}{c}\Delta .Math..Math.X_{\mathrm{dY}}\\ \Delta .Math..Math.P_{\mathrm{dY}}^S\end{array}\right]=-\mathrm{dY};.Math.\left[J_{\mathrm{LF}}^S\right]\left[\begin{array}{c}\Delta .Math..Math.X_{\Delta .Math..Math.F}\\ \Delta .Math..Math.P_{\Delta .Math..Math.F}^S\end{array}\right]=-\Delta .Math..Math.F,$

[0072] where [J.sub.LF.sup.S]=[J, α] is the load flow Jacobian with the distributed slack bus.

[0073] The Lagrange multiplier vector is obtained by solving the following linear equations set:

[00005]${\left[J_{\mathrm{LF}}^S\right]}^T.Math.\mathrm{\Delta \lambda }=-\left[\begin{array}{c}{\nabla }_X.Math.L+\left[H\right].Math.\Delta .Math..Math.X+\mathrm{\beta \Delta }.Math..Math.X\\ 0\end{array}\right].$

[0074] The slack bus power change is determined by the following expression:

ΔP.sup.S=ΔP.sub.dY.sup.SΔt+ΔP.sub.ΔF.sup.S.

[0075] Like the loading factor, the slack bus power is linearly included in (19). For this reason after each iteration it can be corrected by −ΔF.sup.Tα/∥α∥.sub.2.sup.2, thereby compensating the projection of power mismatches onto vector α.

[0076] Let's consider the taking into account a limit-induced bifurcation.

[0077] The treatment of the generator reactive power constraints are performed just as in ordinary load flow solution. However if an MS corresponds to the limit-induced bifurcation (LIB) the iterative process can no converge because of LIB peculiarities: at least one of generator buses switches from PV type to PQ type and vice versa, repeatedly. Such situation occurs when a single generator in achieves its maximum reactive power limit Q.sup.+.sub.m; dQ.sub.m/dt>0, where Q.sub.m is the reactive output of the generator; the Jacobian is not singular, and if generator in switches from PV type to PQ.sup.+ type, the Jacobian determinant changes its sign (when load powers are positive in sign in load flow equations), and dV.sub.m/dt>0. Taking into account such peculiarities allows elaborating the following procedure of detection and determination of LIB.

[0078] The LIB monitoring is performed if the Jacobian determinant changes its sign. Using ΔX from (11), the PQ.sup.+-type generators for which the voltage magnitude at the iteration will exceed the specified value are detected. If such generators are not found, there is no LIB and the iteration of the MS determination procedure continues. On the other hand these generators are switched from PQ.sup.+ to PV type, and the following procedure of LIB determination is performed. [0079] a. Step 1. Solve (12)-(12a) and calculate

[00006]$t=\underset{k.Math.\mathrm{PV},.Math.{\nabla }_X.Math.Q_k^T.Math.\Delta .Math..Math.X_{\mathrm{dY}}>0}{\min }.Math.\left(-\frac{\Delta .Math..Math.Q_k+{\nabla }_X.Math.Q_k^T.Math.\Delta .Math..Math.X_{\Delta .Math..Math.F}}{{\nabla }_X.Math.Q_k^T.Math.\Delta .Math..Math.X_{\mathrm{dY}}}\right),$ [0080] b. where ∇.sub.XQ.sub.k is the gradient of reactive power balance equation for bus k with respect to X, ΔQ.sub.k=Q.sub.k−Q.sup.+.sub.k. [0081] c. Obtain ΔX=ΔX.sub.ΔF+tΔX.sub.dY and update X. [0082] d. Step 2 Check up the constraints for PQ-type generator buses. If constraints violations are found, then switch types of corresponding generators and return to step 1; else continue. [0083] e. Step 3. Compute ΔQ.sub.m=max|ΔQ.sub.k|, k⊂PV. If max{∥ΔF∥.sub.∞, |ΔQ.sub.m|}>ε, then return to step 1; else continue. [0084] f. Step 4. Switch the type of generator bus m from PV to PQ.sup.+ and solve (12). If ΔV.sub.m<0, then return to step 1; else STOP. The LIB is determined: the reactive power limit of generator bus m induces the bifurcation.

[0085] Thus if to take into account the procedure of detection and determination of LIB, the algorithm of the MS determination procedure can be presented as [0086] a. Step 0. Solve the load flow for a base case and obtain an initial guess for λ by the inverse power method. [0087] b. Step 1. Solve (12)-(12a). Obtain Δt.sub.max using ΔX.sub.dY, then Δt and ΔX using (16), (14a), and (11). [0088] c. Step 2. If the Jacobian determinant changes its sign, then run the LIB detection. If the LIB is detected, then run the LIB determination and STOP; else continue. [0089] d. Step 3. Solve (17) to obtain Δλ. Update t, X, and λ, and next correct t by (18). Check for the reactive power limits and switch bus types if needed. [0090] e. Step 4. Convergence check: if ∥ΔF∥.sub.∞<ε and |σ|=∥∇.sub.XL∥.sub.2/∥λ∥.sub.2<ε.sub.σ, then STOP, else return to step 1. [0091] f. The specified algorithm is realized in the proposed device by the following way. [0092] g. In the computerized device for the marginal state assessment of power systems the current network parameters are entered into the operative memory blocks 1-1 . . . 1-1n of the group (FIG. 1) and are collected in a data acquisition block 2 as initial data necessary for performing calculations in a marginal state assessment block 3. The resulted marginal state assessments are put down into an output memory block 4.

[0093] The initial approximation of Lagrange multipliers vector λ is defined for example by the inverse power method in the calculator 5, and the input of this calculator is the input of a marginal state assessment block 3. It allows in a calculator 6 of loading factor ultimate increment by solving the system (12)-(12a) to calculate Δt.sub.max using ΔX.sub.dY and then to calculate Δt ΔX using (16), (14a) and (11). In calculator 7 of power flow matrix determinant sign a check is performed whether a sign of Jacobian matrix determinant has been changed or not. If Jacobian matrix determinant sign has been changed, in a calculator 10 of bifurcation a special procedure of LIB detecting and defining is performed and this fact corresponds to the end of the basic procedure of MS assessment and after that a signal from a calculator 10 output is entered into an output memory block 4. Otherwise a procedure of MS assessment is continued in a calculator 8 of increments and corrections, in this calculator an equation (17) is solved for obtaining Δλ, the values t, X, and λ are modified, and then t is corrected according to (18). Here also the limits on reactive power of generators are checked and if necessary a type of them is changed. After that a signal from the output of calculator 8 is entered to the input of a calculator 9 for a convergence check: ec ∥ΔF∥.sub.∞<ε and σ=∥∇.sub.XL∥.sub.2/∥λ∥.sub.2<ε.sub.σ, then a corrected value of t is entered into an output memory block 4 and it means the end of the basic procedure of MS assessment, otherwise a signal is transferred to a calculator 8 to continue calculations taking into account the modifications and corrections of t, X, and λ.

[0094] Let us consider an example of the proposed device application for two-bus power system. In Table (FIG. 2) the results of its application for MS assessment at the given direction of power change (increase of generation) at bus k are presented. The node data and results are presented for the following parameters of the system: Y.sub.km=−(20+j40).sup.−1Ω.sup.−1, V.sub.k=V.sub.m=110 kV, P.sub.k=100 MW, dP.sub.k=10 MW, α.sub.km=26.565.sup.0, δ.sup.0.sub.km=22.115.sup.0, dΔP.sub.k/dδ.sup.0.sub.km−269.748.

[0095] The system has two MS. One of them corresponds to the maximum generation at bus k, another—to the maximum load at this bus. The maximum generation turns out when δ.sub.km−α.sub.km=π/2, thus the injected power will be equal to P.sub.k.sup.max gen=V.sub.k.sup.2|Y.sub.km|sin α.sub.km+V.sub.kV.sub.m|Y.sub.km|=391.567 MW. The maximum load turns out when δ.sub.km−α.sub.km=−π/2, thus the consumed power will be equal to P.sub.k.sup.max load=V.sub.k.sup.2|Y.sub.km|sin α.sub.km−V.sub.kV.sub.m|Y.sub.km|=−149.567 MW at bus k. In the border-state (BS) i.e., when δ.sub.km=α.sub.km, bus k injects P.sub.k.sup.BS=V.sub.k.sup.2|Y.sub.km|sin α.sub.km=121 MW. Thus, the BS appears to be the equidistance state from these two MS.

[0096] Bus k injects 100 MW in the base case. Since P.sub.k<P.sub.k.sup.BS, this state is close to the second MS. In the base case δ.sub.km.sup.0=22.115.sup.0 and the line loss angle α.sub.km=26.565.sup.0. Since δ.sub.km<α.sub.km, therefore the system solution (3)-(3b) should be the MS in the opposite direction of the given power change. I.e., model (1)-(1a) moves the system state to the closest MS in the power space. Moreover, the system matrix (4) is ill-conditioned in the initial state. The value of denominator (14) is equal to ΔX.sub.dY.sup.T[H]ΔX.sub.dY=−2.885.Math.10.sup.−3. If β is not used, then the system solution (4) will be Δt.sup.1=−346.595 and Δδ.sub.km.sup.1=−736.1830.sup.0=−(4π+16.183.sup.0) at the first iteration. I.e., the voltage angle δ.sub.km will be turned more than two times in the “wrong” direction. The proposed system overcomes these difficulties easily. FIG. 2 shows that ΔX.sub.dY.sup.T[H]ΔX.sub.dY is negative at the first iteration. Therefore, according to (16) the matrix [H] is “corrected” by means of β. In order to improve the condition number, at each iteration the value of Δt.sub.max is determined so that the change of the voltage angle (ΔX.sub.dYΔt) across the line does not exceed 45.sup.0. At the second iteration the value of ΔX.sub.dY.sup.T[H]ΔX.sub.dY becomes positive, and β is used for improvement of the condition number. The MS calculation requires 3 iterations only. After the third iteration the power mismatch |ΔP.sub.k| is equal to 3.189.Math.10.sup.−3 MW, and the load flow Jacobian eigenvalue achieves almost zero value, namely 1.032.Math.10.sup.−5. Thus, the proposed system provides the fast, reliable and exact MS calculation. The accuracy is confirmed easily by the expression δ.sub.km.sup.MS=π/2+α.sub.km=90.sup.0+26.565.sup.0=116.565.sup.0 which should be carried out in the MS. After 3 iterations the voltage angle δ.sub.km is equal exactly to this value.

[0097] Thus, due to the introduction of an additional arsenal of technical facilities the required technical result consisting in an increase of the device speed of response when performing the MS assessment in power systems, is achieved that is confirmed by experiment. The additional arsenal of technical facilities particularly consists in the following: the outputs of a group of operative memory blocks are connected with the input of a data acquisition block and the output of this block is connected with the input of a marginal state assessment block with its output connected to the input of an output memory block, a marginal state assessment block is performed in a form of connected in series a calculator of Lagrange multiplier vector with its input being the input of a marginal state assessment block, a calculator of loading factor ultimate increment, a calculator of a sign of power flow matrix determinant, a calculator of increments and corrections and a calculator of convergence check that has the first output which is the first output of a marginal states assessment block and the second output connected to the additional input of a calculator of loading factor ultimate increment, and also a calculator of bifurcation with the input connected to the additional output of a calculator of a sign of power flow matrix determinant and the output which is the second output of a marginal state assessment block, besides the first and the second outputs of a marginal state assessment block are connected correspondingly with the first and the second inputs of an output memory block.