Space Division and WGAN-GP Based Fast Generation Method of Practical Dynamic Security Region Boundary

Yanli Liu; Ruipeng Jia

doi:10.1016/j.eng.2024.05.017

Engineering ›› 2025, Vol. 51 ›› Issue (8) :80 -91. DOI: 10.1016/j.eng.2024.05.017

Research

research-article

Space Division and WGAN-GP Based Fast Generation Method of Practical Dynamic Security Region Boundary

Yanli Liu ^*
, Ruipeng Jia

Author information +

History +

PDF (4940KB)

Abstract

Fast and accurate transient stability analysis is crucial to power system operation. With high penetration level of wind power resources, practical dynamic security region (PDSR) with hyper plane expression has outstanding advantages in situational awareness and series of optimization problems. The precondition of obtaining accurate PDSR boundary is to locate sufficient points around the boundary (critical points). Therefore, this paper proposes a space division and Wasserstein generative adversarial network with gradient penalty (WGAN-GP) based fast generation method of PDSR boundary. First, the typical differential topological characterizations of dynamic security region (DSR) provide strong theoretical foundation that the interior of DSR is hole-free and the boundaries of DSR are tight and knot-free. Then, the space division method is proposed to calculate critical operation area where the PDSR boundary is located, tremendously compressing the search space to locate critical points and improving the confidence level of boundary fitting result. Furthermore, the WGAN-GP model is utilized to fast obtain large number of critical points based on learning the data distribution of the small training set aforementioned. Finally, the PDSR boundary with hyperplanes is fitted by the least square method. The case study is tested on the Institute of Electrical and Electronics Engineers (IEEE) 39-bus system and the results verify the accuracy and efficiency of the proposed method.

Graphical abstract

Keywords

Practical dynamic security region boundary / WGAN-GP / Space division / Critical points / Data-model hybrid driven

Cite this article

Download citation ▾

Yanli Liu, Ruipeng Jia. Space Division and WGAN-GP Based Fast Generation Method of Practical Dynamic Security Region Boundary. Engineering, 2025, 51(8): 80-91 DOI:10.1016/j.eng.2024.05.017

登录浏览全文

4963

注册一个新账户忘记密码

1. Introduction

The access of high percentage of renewable energy sources brings great uncertainty to the power system [1], [2] and leverage the power system stability facing massive risks and challenges [3], [4]. Recently, large-scale power blackout accidents have been closely associated with power system transient instability [5]. Therefore, fast and accurate transient stability [6], [7] analysis is urgently needed to ensure secure operation of power systems.

Point-by-point method [8], which conducts time-domain simulation to judge stability results point by point, has been long-time used to analyze the security and stability of power system. Dynamic security region (DSR), developed on the basis of point-by-point method, is defined in the node injection space before the accident, and refers to the set of points which can still satisfy the transient power angle stability of the system after the accident. DSR implements offline calculation and online application, which not only judges the system stability, but obtains the security margin and control direction based on the relative position between the operating point and the DSR boundary [9]. The DSR boundary can be approximated by one or a few hyperplanes (HPs), called practical DSR (PDSR) [10]. Owing to the availability of convenient and concise expressions, the PDSR is a powerful tool for solving optimization problems with transient power angle stability constraints and power injection uncertainties, providing schedulers with plentiful system state information and control decision reference [11], [12] and having good prospects for engineering applications.

Traditionally, the typical methods for the PDSR boundary generation include analytical methods [13] and fitting methods [14]. The analytical method derives an analytical formula of the PDSR boundary based on the mapping relationship between the node power injection and system transient stability state. Ref. [13] derives a PDSR boundary analytical expression which is more convenient to calculate based on the property that the security region boundaries are parallel to each other with small injection changes. Ref. [15] quickly identifies the instability modes of the critical points around the initial operating point and provides suggestions for critical points searching of the analytical method. The fitting method fits the expression of PDSR boundary based on a large number of critical points by time-domain simulation, where critical points are the operating points of the transition state between the stable and unstable states. Ref. [14] obtains the PDSR boundary based on the orthogonal point selection method and the least square method, which lays the foundation for many relevant researches. Ref. [16] applies PDSR to the actual large scale power system of Central China Branch of State Grid Corporation of China, proving the practicality of PDSR boundary and its relevant properties. The above researches show that the analytical methods only search very few critical points with fast calculation speed, but fail to meet high calculation accuracy because of the power system complex model. Whereas, the fitting methods focus on calculation accuracy performance, which is done at the expense of time consumption (i.e., with increasing scale of the power system, higher power injection space dimensions require more work on time-domain simulation). Therefore, how to obtain PDSR boundary considering accuracy and efficiency simultaneously plays a crucial role in practical engineering application.

Data driven or hybrid data driven method for obtaining the boundary is proposed recently. Relief method [17] is utilized to identify the key generators and reduce the data dimension. Generative adversarial network (GAN) [18], [19], [20], [21] is used for date augment and to obtain boundaries based on the distribution of massive stable and unstable points. Refs. [22], [23], [24] apply advanced deep learning methods to evaluate transient stability, which rely on time-domain simulation to generate numerous operating points with corresponding transient stability labels as train sets. However, it is time consuming to generate the datasets. Since the operating points farther from the boundary contain less information about transient stability state, it is more efficient to generate operating points close to the boundary rather than within the whole injection space. Nevertheless, it lacks relevant demonstration about the confidence region of critical point distribution, that is there is no guarantee that generated points are exactly distributed on the boundary or very close to it.

This paper proposes a space division and Wasserstein GAN with gradient penalty (WGAN-GP) based method for PDSR boundary fast generation. First, the typical differential topological characterizations of DSR are denseness, continuity, and compactness, which provide strong theoretical foundation for searching for critical points only in a small space instead of the whole injection space. The space division method is promoted to compress the search space by dividing the injection space into three areas. A small set of critical points can be searched only in the critical operating area where the PDSR boundary is located. Then, with Wasserstein distance and gradient penalty as the loss function, the WGAN-GP model is utilized to generate large number of critical points in the critical operating area, realizing data augment of the critical points aforementioned. Finally, the PDSR boundary with HPs is fast fitted based on the least square method, which enables efficient and accurate power system transient stability assessment.

The remainder of this paper is organized as follows. Section 2 formulates the definition of DSR and PDSR and provides the typical differential topological characterizations of DSR. Section 3 introduces the space division method. Section 4 addresses the implementation of the proposed method. The experimental results and discussions are shown in Section 5, and conclusions are drawn in Section 6.

2. Dynamic security region

2.1. Dynamic security region

DSR is the set of points in the pre-accident injection space, which can guarantee the transient stability of the power system after the accident. In high-voltage alternating current (AC) power grid, it can be assumed that the reactive power can be balanced in location, and the effect of the change of active power on the magnitude of the voltage is very small and negligible. Then, the DSR can be studied only in the space of the active power injection. The change in state of the system from the occurrence of a disturbance to its removal can be represented by a set of differential equations.

(1)

{\dot{x}}_{β} (t) = \{\begin{matrix} f_{i} (x_{β} (t)), \\ f_{F} (x_{β} (t)), \\ f_{j} (x_{β} (t)), \end{matrix} \begin{matrix} - \infty < t < 0 \\ 0 \leq t < τ \\ τ \leq t < + \infty \end{matrix}

where i and j are the power system structure before and after the disturbance, respectively; F denotes the power system structure in the disturbance; t denotes the time phase of the power system;

x_{β} (t)

denotes injection power vector at time t;

τ

is the duration of the disturbance; and f(·) denotes the mapping relation from the decision space to the state of the system;

{\dot{x}}_{β} (t)

denotes the transient stable state of the power system at time t. Specifically, the DSR

Ω_{d} (i, j, τ)

is defined as follows.

(2)

Ω_{d} (i, j, τ) ≜ {x_{β} | x_{d} (x_{β}) \in A (x_{β})}

where

x_{d} (x_{β})

is the system state at the moment of disturbance clearance;

A (x_{β})

is the transient stability region of the equilibrium point determined by the injection power vector

x_{β}

in the state space after the disturbance. The mapping relationship between the DSR

Ω_{d} (i, j, τ)

and the transient stability region

A (x_{β})

is shown in Fig. 1.

2.2. Typical differential topological characterizations of DSR

For the power system structure preserving model [25], the system state can be expressed as

x = (δ, ω)

, where

δ

denotes the phase angle vector of the load bus and generator bus of the system, and

ω = \dot{δ}

denotes the frequency deviation vector of the bus. The gradient system is defined as follows.

(3)

\dot{δ} = - \frac{∂ V_{p} (δ)}{∂ δ}

where

V_{p}

represents the potential energy function of the system.

According to Eq. (2), define the DSR

Ω_{g} (i, j, τ)

of the gradient system.

(4)

Ω_{g} (i, j, τ) ≜ {x_{β} | x_{g} (x_{β}) \in A_{g} (x_{β})}

(5)

x_{d} (x_{β}) = (x_{g} (x_{β}), ω)

where

x_{g} (x_{β})

is obtained by the projection of the system instantaneous state

x_{d} (x_{β})

at the time of the accident clearance on the

δ

-space.

A_{g}

indicates the stability region of the original system and

∂ A_{g}

denotes the stability region boundary which is also known as the potential energy surface.

A_{g} (x_{β})

is the stability region corresponding to the gradient system.

The nonlinear autonomous dynamic system satisfies the following assumptions: ① All equilibrium points on the stability boundary are hyperbolic; ② both stable and unstable manifolds of the equilibrium point on the stability boundary satisfy the transversality condition; ③ the number of equilibrium points on the stability boundary is finite; ④ every trajectory of the stability boundary tends to an equilibrium point as time goes to infinity.

Then, it can be deduced that the gradient system of power system is structurally stable [26]. Based on the special properties of n-dimensional gradient systems and the potential energy boundary surface (PEBS) approach, the differential topological characterizations of the DSR of original system can be proved by the DSR of gradient system.

(1) For any point y belonging to

Ω_{g} (i, j, τ)

, there is always a unique gradient system

\dot{δ} = - ∂ V_{p} (δ) / ∂ δ

corresponding to the point, which is first-order continuous differentiable, expressed as C¹. Thus, there exists a mapping

β (\cdot) : x_{β} \to - ∂ V_{p} / ∂ δ, \forall x_{β} \in Ω_{g} (i, j, τ)

, which is continuous and

β^{- 1} (\cdot)

is also continuous. Thus,

β (\cdot)

is homeomorphism. In this case, the set of gradient systems and the corresponding DSR

Ω_{g} (i, j, τ)

have the same topology under homeomorphism

β (\cdot)

Thus, it can be deduced that the gradient systems corresponding to

x_{β}

belongs to

Ω_{g} (i, j, τ)

, are structurally stable.

Moreover, at any point within the boundary of the corresponding stability region of

x_{β}

, the gradient vector field

f (δ) = - ∂ V_{p} (δ) / ∂ δ

points inward along the boundary of the stability region. According to the structural stability theory, the set of gradient systems corresponding to

x_{β}

is dense and the DSR

Ω_{g} (i, j, τ)

is also dense.

In order to further demonstrate that the DSR

Ω_{d} (i, j, τ)

of the original system is dense, an approximate description of the local boundary near the kth unstable equilibrium point (UEP)

({\hat{δ}}_{k}, 0)

on the stability boundary of the original system can be used to obtain the approximate DSR, which is defined as follows.

(6)

{\tilde{Ω}}_{d} (i, j, τ) ≜ {x_{β} | x_{d} (x_{β}) \in \tilde{A} (x_{β})}

where,

\tilde{A} (x_{β})

is a stability region determined by the energy of the relevant UEP as the critical energy, expressed as follows.

(7)

∂ \bar{A} (x_{β}) ≜ ⋃_{k} ∂ \bar{V} ({\hat{δ}}_{k}, 0)

where

∂ \bar{V} ({\hat{δ}}_{k}, 0) ≜ \{(δ, ω) | V (δ, ω) = V ({\hat{δ}}_{k}, 0) = V_{p} ({\hat{δ}}_{k})\}

. Then, proving that the approximate DSR is dense can verify the DSR of the original system is dense.

Assume that the trajectory of the accident tends to the stable manifold of the kth UEP

({\hat{δ}}_{k}, 0)

on the boundary of the original system. According to Eqs. (6), (7) and relevant theory of the energy function,

\forall x_{β 1} \in {\tilde{Ω}}_{d} (i, j, τ)

, Eq. (8) is established.

(8)

V_{p} (x_{g} (x_{β 1})) = V (x_{g} (x_{β 1}), 0) < V (x_{d} (x_{β 1})) < V ({\hat{δ}}_{k}, 0) = V_{p} ({\hat{δ}}_{k})

Since the potential energy function

V_{p}

of the original system is the energy function of the gradient system, the energy

V_{p} ({\hat{δ}}_{k})

at the relevant UEP is the smallest on the stable boundary of the gradient system. Therefore,

x_{g} (x_{β 1})

belongs to

A_{g}

(i.e.,

x_{β 1}

belongs to

Ω_{d} (i, j, τ)

). And

{\tilde{Ω}}_{d} (i, j, τ) \subset Ω_{g} (i, j, τ)

{\tilde{Ω}}_{d} (i, j, τ)

is also dense.

In addition, it is considered that the method based on UEP in transient stability analysis of power system usually has an acceptable accuracy range. The approximate DSR determined by the correlation UEP method is close to the real DSR, so it can be inferred that the DSR

Ω_{d} (i, j, τ)

is open and dense.

(2) The essence of the phenomenon of “knotting” on the boundary of DSR is that the manifold on the boundary has a limit cycle or chaotic motion. However, since the boundary of the stability region satisfies the assumptions ①–④, it can be seen that the system will not have bifurcation. There will be no bounded oscillation like the limit cycle or chaotic motion. Therefore, the boundary of DSR is knot-free and the DSR

Ω_{d} (i, j, τ)

is connected.

(3) According to (1), the DSR is open, and the boundary of the open set is a closed set. Thus, the boundary of DSR is a closed set. Since the injection power is within the upper and lower constraints, the boundary is bounded. Therefore, the boundary of DSR is compact.

It can be concluded that the typical differential topological characterizations of DSR are denseness, continuity, and compactness, that is the interior of DSR is hole-free and the boundaries of DSR are tight and knot-free. And the proof detail can be referred to the Refs. [27], [28]. The typical characterizations lay a solid theoretical foundation for searching for critical points only in a small space instead of the whole injection space.

2.3. Practical dynamic security region

In practical engineering, it is assumed that the reactive power of the high-voltage AC power system is balanced locally. Then, given system topology, accident type, and clearance time, DSR can be approximated by one or a few HPs, called PDSR. In n-dimensional space, the above HPs can be depicted as follows.

(9)

\sum α_{v} P_{v} = 1, \forall v \in (1, . . ., n)

where

(α_{1}, . . ., α_{n})

is the coefficient vector;

(P_{1}, . . ., P_{n})

is the active injection power vector of critical points. According to the characteristics of PDSR, it is generally considered that the operating point satisfying

\sum α_{v} P_{v} < 1

is transiently stable, while the one satisfying

\sum α_{v} P_{v} > 1

is transiently unstable.

3. Space division method

Provided the differential topological characterizations of DSR with denseness, continuity, and compactness, space division method precisely determines the critical operating area where critical points are concentrated in, owing to the fact that critical points are not uniformly distributed in the whole injection space, but in a specific area. Then, a large number of critical points can be fast obtained just in the critical operating area instead of the whole injection space, which is time-saving by compressing the search space.

Space division method divides the node active power injection space into three parts, including stable operating area with concentrated distribution of stable operating points, unstable operating area with concentrated distribution of unstable operating points, and critical operating area with concentrated distribution of critical operating points and the PDSR boundary. The schematic of the divided node injection space is shown in Fig. 2, where S₀ indicates the reference point, crP is the set of critical points, and D_l represents the average Euclidean distance between the reference point and the set of critical points (l).

[a \times D_{l}, b \times D_{l}]

denotes the critical operating area, where a and b are coefficients which guarantee the engineering error is within a reasonable range.

Since the critical points are within the critical operating area, the initial critical point determined by dichotomy and time-domain simulation can be utilized to calculate the critical operating area, guiding the acquisition of more critical points. Moreover, the reference distance from the reference point to the boundary is roughly equivalent to the distance between the reference point and the critical points. Then, the critical operating area can be determined by computing the distance between the reference point and the critical points. The flowchart for calculating critical operating area is shown in Fig. 3 and the specific calculation method is as follows.

The Euclidean distance d_q between the reference point S₀ and the qth critical point crP_q is calculated as follows.

(10)

d_{q} = \sqrt{{(S_{0_G 1} - cr P_{q_G 1})}^{2} + . . . + {(S_{0_L 1} - cr P_{q_L 1})}^{2} + . . .}

where G1 is the first generator and L1 is the first line load.

Take the Euclidean distance d₁ between the reference point S₀ and the first critical point crP₁ as the initial reference distance D₁, and the critical operating area can be derived as

[a \times D_{1}, b \times D_{1}]

, where a is set to 0.8 and b is set to 1.2.

When a new critical point is searched, the critical operating area

[a \times D_{l}, b \times D_{l}]

is updated based on the average Euclidean distance D_l between the reference point and l critical points that already searched.

(11)

D_{l} = (d_{1} + . . . + d_{n}) / n

Calculate the Euclidean distance d₀ between the reference point S₀ and the operating point P.

(12)

d_{0} = \sqrt{{(S_{0_G 1} - P_{G 1})}^{2} + . . . + {(S_{0_L 1} - P_{L1})}^{2} + . . .}

During the application process, by comparing the values of d₀ and D_l, the location relationship between the operating point and the critical operating area can be determined. Then, space division method can guide the searching for critical points in the critical operating area. If the operating point is located within the critical operating area, continue searching and calculating based on the dichotomy and time-domain simulation. If the operating point is located outside the critical operating area, it indicates that the operating point is not a critical point. Then change the injection power and calculate the next operating point.

4. Proposed methodology

The focus of PDSR boundary fast generation is to search sufficient critical points among massive operating points. Meanwhile, the precision of the critical points guarantees the high-accuracy of the fitted boundaries. However, the expansion of the system scale explodes the number of critical points to be obtained, continuously increasing the calculation burden.

To address the limitation, the framework of the proposed method is illustrated in Fig. 4. Under given scenario, accident type, and clearance time, calculate the first critical point and obtain the initial operating area. Then space division method is promoted to compress search space for calculating a small training set of critical points just in the critical operating area. Additionally, WGAN-GP model is utilized to fast-generate numerous critical points in the critical operating area by learning the data distribution of the set of critical points. Finally, the PDSR boundary with HPs is fitted by the least square method. The key steps are elaborated in the sequel.

4.1. Fast search for a small number of critical points

A small number of critical points are searched by space division and Relief to compress the search space and reduce the feature dimension. And to ensure reliable samples and accurate calculations, obtain at least one-third of the required number of critical points. The flowchart of critical points searching method is shown in Fig. 5 and the key steps is as follows.

Step 1: Identify the key units closely related to the system stability based on the Relief algorithm and the massive operating points under the given fault. Due to the page limit, the detail can be referred to the Ref. [29];

Step 2: Determine the initial search space of the critical points according to the Hadamard matrix, and take each row of the orthogonal matrix as the search operating point P. According to the dichotomous method and time-domain simulation, search and calculate until the initial critical point crP₁ is obtained. Based on the space division method, compute the initial critical operating area;

Step 3: Determine the extended search space based on the initial critical point crP₁ and the Hadamard matrix, and take each row of the orthogonal matrix as the search operating point P. If P is in the critical operating area

[a \times D_{l}, b \times D_{l}]

, search the critical point based on the dichotomy and time domain simulation, where

D_{l}

is the reference distance when the lth critical point is searched. If P is outside the area, change the operating point and continue searching;

Step 4: Update the critical operating area based on the average Euclidean distance

D_{l}

between the reference point S₀ and the searched critical points. Return to step 3 until s critical points are searched.

4.2. Critical point augment based on WGAN-GP

Sufficient critical points can be fast generated by data augment based on a small number of critical points. GAN is an efficient network for data augment and a variant of GAN which can overcome the problems of gradient disappearance and model collapse of the original GAN is urgently needed.

By introducing Wasserstein distance and gradient penalty, WGAN-GP offers unique advantages in enhancing training stability and improving model effect with gradient disappearance and model collapse problems solved. Therefore, WGAN-GP model is adopted to fast generate a large number of critical points based on the injected power of critical points obtained in Section 4.1.

The Wasserstein distance in WGAN-GP model is defined as follows.

(13)

W (P_{data}, P_{G}) = \inf_{γ \sim Π_{(P_{data}, P_{G})}} E_{(x, y) \sim γ} [‖ x - y ‖]

where P_data is the true sample distribution (i.e., the critical point distribution); P_G is the generated sample distribution;

Π_{(P_{d a t a}, P_{G})}

is the set of P_data and P_G joint distributions;

γ

is one of the joint distributions

Π_{(P_{d a t a}, P_{G})}

;

E_{(x, y) \sim γ}

denotes the expectation, where (x, y) obeying the distribution

γ

‖ x - y ‖

is the norm of

x - y

;

W (P_{data}, P_{G})

is the lower bound of the expectation of

‖ x - y ‖

among all the possible distributions.

Since

\inf_{γ \sim Π_{(P_{data}, P_{G})}}

cannot be solved directly, the Lipschitz continuum is introduced to convert Wasserstein distance into Eq. (14).

(14)

W (P_{data}, P_{G}) = \frac{1}{K} \sup_{{‖ f ‖}_{L} \leq K} E_{x \sim P_{data}} [f (x)] - E_{x \sim P_{G}} [f (x)]

where

f (x)

is the equivalent mapping function of the discriminator;

{‖ f (x) ‖}_{L} < K

denotes that the Lipschitz constant of

f (x)

does not exceed K (i.e.,

f (x)

satisfies Lipschitz continuity); K is a constant which is not positive infinity.

The objective function of WGAN-GP model is shown as follows.

(15)

L = \max E_{x \sim P_{data}} [D_{μ} (x)] - E_{x \sim P_{G}} [D_{μ} (x)] - λ \times GP

(16)

GP = E_{\hat{x} \sim p (\hat{x})} \max (∥\nabla_{\hat{x}} D_{μ} (\hat{x})∥, 1)

where

D_{μ} (x)

is the output of the discriminator

D_{μ}

with upper bound variables

μ

when input x. Point

\hat{x}

satisfies

\hat{x} = ϖ x^{data} + (1 - ϖ) x^{G}

ϖ \sim u (0, 1)

, where

x^{data}

is the real sample,

x^{G}

is the generated sample;

λ

is the penalty factor; GP is the gradient penalty;

∥\nabla_{\hat{x}} D_{μ} (\hat{x})∥

denotes the norm of gradient

\nabla_{\hat{x}} D_{μ} (\hat{x})

Then, the loss functions of the generator and the discriminator can be expressed as Eqs. (17), (18), respectively.

(17)

L_{G} = - E_{x \sim P_{G}} [D_{μ} (x)]

(18)

L_{D} = E_{x \sim P_{data}} [D_{μ} (x)] - E_{x \sim P_{G}} [D_{μ} (x)] + λ E_{\hat{x} \sim p (\hat{x})} \max (∥\nabla_{\hat{x} \sqrt{b^{2} - 4 a c}} D_{μ} (\hat{x})∥, 1)

During the training process, real critical points are normalized and marked as true first, with generated samples marked as false. Then, the generator implements a mapping from random noise distribution to generated sample distribution. And the discriminator receives both types of samples with corresponding labels and outputs the probability of true for a given sample. In adversarial alternating training, based on the above loss functions, the Adam optimizer is used to update the parameters of the generator and discriminator. With updated parameters, the generator can generate samples that are similar to real critical points and indistinguishable by the discriminator, and the discriminator can accurately distinguish whether a given sample is real or not.

During the application process, a group of random noises subject to simple normal distribution are input into the generator, which can generate a large number of samples that follow the distribution of critical point. Therefore, based on a small number of critical points, the WGAN-GP model can provide an accurate and sufficient sample set of critical points for specific accidents.

4.3. Fit the PDSR boundary

Based on the least square method, the PDSR boundary can be obtained by fitting. Assuming that in the n-dimensional node injection space, the required sample set includes k critical points. Then, the active injection vector of the qth critical point is denoted as

P_{q} = (P_{q 1}, . . ., P_{qn})

, and the associated deviation equation can be depicted as Eq. (19).

(19)

Y - P α = ε

(20)

P = [\begin{matrix} P_{11} & P_{12} & \dots & P_{1 n} \\ P_{21} & P_{22} & \dots & P_{2 n} \\ ⋮ & ⋮ & ⋮ \\ P_{k 1} & P_{k 2} & \dots & P_{kn} \end{matrix}]

where

Y = {[1, . . ., 1]}^{T}

is a vector of order

k \times 1

, representing measured value; P is a matrix of order

k \times n

with constant elements;

α = {[α_{1}, . . ., α_{n}]}^{T}

is a vector of order

n \times 1

, representing the HP coefficient corresponding to each node; the deviation

ε = {[ε_{1}, . . ., ε_{k}]}^{T}

is a vector of order

k \times 1

, denoting fitting error.

Based on the theory of the least square method, the sum of squares of deviations should be minimized, as shown in Eq. (21).

(21)

J = \min {∥ε∥}^{2} = \min {(Y - P α)}^{T} (Y - P α)

Then, calculate the extreme value by partial derivative as shown in Eq. (22).

(22)

\frac{∂ J}{∂ α} = \frac{∂ (Y^{T} Y - 2 Y^{T} P α + α^{T} P^{T} P α)}{∂ α} = - 2 Y^{T} P + 2 α^{T} P^{T} P = 0

After solving the above equations, the least square estimate of the HP coefficient can be depicted as shown in Eq. (23).

(23)

α_{Ls} = {(P^{T} P)}^{- 1} P^{T} Y

The distance between the critical point and the edge of the PDSR boundary with HPs is taken as the fitting error

err = {[{err}_{1}, \dots, {err}_{k}]}^{T}

. The fitting error

{err}_{q}

for the qth critical point is calculated as shown in Eq. (24).

(α_{1}, . . ., α_{v})

is the coefficient vector;

(X_{q 1}, . . ., X_{qv})

is the active injection power vector of critical points and v represents the number of features. The maximum fitting error is defined as

{err}_{m} = \max (err, . . ., {err}_{k})

to represent the approximate calculation accuracy of the PDSR boundary. The smaller the maximum fitting error is, the higher the approximate calculation accuracy is, and vice versa.

(24)

{err}_{q} = \frac{|\sum_{i = 1}^{n} α_{i} X_{qi} - 1|}{\sqrt{\sum_{i = 1}^{n} α_{i}^{2}} \cdot \sqrt{\sum_{i = 1}^{n} X_{qi}^{2}}}

5. Case study

5.1. Setting of the case study

The Institute of Electrical and Electronics Engineers (IEEE) 39-bus test system with generator 35 replaced by a wind turbine is chosen to verify the efficiency of the proposed method shown in the left part of Fig. 6. In order to enrich the case study, The IEEE 39-bus test system with generator 35 and generator 38 separately replaced by a wind turbine is also chosen to verify the proposed method shown in the right part of Fig. 6. Generator 31 bus is selected to be the slack bus. Fault line 1–2 denotes a three-phase short-circuit fault at bus 1 of line 1–2. Table 1, Table 2 show the relevant data of the IEEE 39-bus test system. In order to simplify the description, the following analysis is based on the IEEE 39-bus test system with generator 35 replaced by a wind turbine, except for the two wind turbines system mentioned. The proposed method is tested on the computer with 11th Gen Intel(R) Core(TM) i5-11300H central processing unit (CPU)@3.10 GHz with 3.11 GHz mains frequency and random access memory (RAM) 16.00GB. The transient simulation is calculated using Bonneville Power Administration (BPA) software. The algorithmic parts are calculated on MATLAB and python.

5.2. Case analysis

Table 3 illustrates that the points in the stable operating area are of stable states and Table 4 illustrates that the points in the unstable operating area are of unstable states. Moreover, Table 5 strengthens that space division method divides the operating area precisely which leads to high accuracy of transient stability of operating points in the corresponding area.

Table 6 shows that the maximum fitting error of the space division method is on the order of 1.0

\times

10⁻⁴, indicating high calculation accuracy.

Fig. 7 exhibits that the space division method reduces calculation time by 71.74% at most, 17.96% at least, and 34.54% on average compared to the fitting method, emphasizing the advantage of compressing the search space to reduce time consumption.

In order to ensure fitting accuracy, the number of critical points should be greater than the dimension of the PDSR injection space with the quasi-orthogonal point selection method adopted [30]. The injection space of the IEEE 39-bus test system is 28 dimensions, so more than 28 critical points need to be searched to generate the PDSR boundaries in this system. To guarantee calculation accuracy and statistical significance, 840 samples are generated based on the WGAN-GP model, which is 30 times the required critical points.

Table 7 shows that the samples generated by WGAN-GP are in the critical operating area and conform to the distribution rule of critical points. It proves that the generated critical points are accurate and effective, laying a foundation for the accurate generation of the PDSR boundary.

Fig. 8 shows the PDSR boundaries for faults 1–2, where the stars represent the critical points and the blue planes represent the PDSR boundary fitted by the critical points. Below the boundary is the PDSR. Since the injection is a multi-dimensional vector, the DSR can only display a certain cross-section in three-dimensional space, for which the injection power of other dimensions must be given, and the location of the boundary is different for different given values of other dimensions.

Table 8 shows that the maximum fitting error of the proposed method is on the order of 1.0 × 10⁻⁵, and less than that of fitting method. Table 9 shows that the maximum fitting error of the proposed method is on the order of 1.0 × 10⁻⁴, and less than that of fitting method in the two wind turbines system, demonstrating high calculation precision.

Fig. 9 exhibits that the long execution time of the proposed method is reduced by 78.76% at most, 40.41% at least, and 55.01% on average compared with the space division and Relief method which has already improved the calculation speed. Table 10 exhibits that the long execution time of the proposed method is reduced by 72.16% on average compared with the fitting method and 61.73% on average compared with the space division and Relief method in the two wind turbines system. Moreover, it is noteworthy that the long execution time of PDSR boundary generation by the proposed method is less than 6 min. Therefore, the proposed method reduces calculation time visually and significantly.

6. Conclusions

This paper proposes a PDSR boundary generation method based on space division and WGAN-GP for fast obtaining sufficient critical points. Provided the typical differential topological characterizations of DSR, space division is promoted to compress the search space and search for a small set of critical points in the critical operating area. Leveraging WGAN-GP model to learn the data distribution of the set aforementioned, plentiful critical points can be fast generated just by one training, and then the PDSR boundary can be fitted precisely. The case study is conducted in the IEEE 39-bus test system with generator 35 replaced by a wind turbine and the IEEE 39-bus test system with generator 35 and generator 38 replaced separately by a wind turbine to validate the efficiency of the proposed method.

Test results show that comparing with the fitting method and its improved algorithm, the speed of the proposed method to obtain numerous critical points and calculate PDSR boundaries is improved dramatically, while allowing for high calculation accuracy.

The proposed method is aimed at given faults. Considering that there is a certain correlation between different faults, with the development of deep learning, it will be the focus of future research to generate PDSR boundaries corresponding to other faults from given faults.

CRediT authorship contribution statement

Yanli Liu: Writing – review & editing, Visualization, Supervision, Methodology, Formal analysis. Ruipeng Jia: Writing – original draft, Visualization, Methodology.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was funded in part by the National Key Research and Development Program of China (2020YFB0905900), and in part by Science and Technology Project of State Grid Corporation of China (SGCC): The Key Technologies for Electric Internet of Things (SGTJDK00DWJS2100223).

References

Publishing order | Descend order by publishing year | Descend order by cited within

[1]	Holttinen H, Kiviluoma J, Flynn D, Smith JC, Orths A, Eriksen PB, et al.System impact studies for near 100% renewable energy systems dominated by inverter based variable generation.IEEE Trans Power Syst 2022; 37(4):3249-3258.

[2]	Liang X.Emerging power quality challenges due to integration of renewable energy sources.IEEE Trans Ind Appl 2017; 53(2):855-866.

[3]	Papadopoulos PN, Milanovic JV.Probabilistic framework for transient stability assessment of power systems with high penetration of renewable generation.IEEE Trans Power Syst 2017; 32(4):3078-3088.

[4]	Xia S, Ding Z, Shahidehpour M, Chan KW, Bu S, Li G.Transient stability-constrained optimal power flow calculation with extremely unstable conditions using energy sensitivity method.IEEE Trans Power Syst 2021; 36(1):355-365.

[5]	Liang Z, Ge R, Dong Y, Chen G.Analysis of large-scale blackout occurred on July 30 and July 31, 2012 in India and its lessons to China’s power grid dispatch and operation.Power Syst Technol 2013; 37(7):1841-1848.

[6]	Hatziargyriou N, Milanovic JV, Rahmann C, Ajjarapu V, Canizares C, Erlich I, et al.Definition and classification of power system stability—revisited & extended.IEEE Trans Power Syst 2021; 36(4):3271-3281.

[7]	Yang T, Yu Y.Static voltage security region-based coordinated voltage control in smart distribution grids.IEEE Trans Smart Grid 2018; 9(6):5494-5502.

[8]	Yu Y.Methodology of security region and practical results.J Tianjin Univ 2003; 36(5):525-528.

[9]	Wang X, Zhang M, Deng M, Qiao Y, Hu X.Optimal emergency control strategy algorithm of ideal time based on dynamic security region.Power System Protect Control 2014; 42(12):71-77.

[10]	Feng F, Yu Y.Dynamic security regions of power systems in injection spaces.Proc CSEE 1993; 13(3):16-24.

[11]	Yu Y, Liu H, Zeng Y.Optimal transient stability emergency control based on practical dynamic security region.Sciencepaper Online 2004; 34(5):556-563.

[12]	Liu H, Yu Y.A comprehensive security control method based on practical dynamic security regions of power systems.Proc CSEE 2005; 25(20):31-36.

[13]	Yu Y, Lin J.Practical analytic expression of power system dynamic security region’s boundary.J Tianjin Univ 1997; 30(1):2-9.

[14]	Yu Y, Luan W.Determination of dynamic security region of practical power system by fitting technology.Proc CSEE 1990; 10(S1):24-30.

[15]	Liang M, Yu Y, Stephen TL, Pei Z.Identification method of instability modes and its application in dynamic security region.Autom Electr Power Syst 2004; 28(11):28-32.

[16]	Zeng Y, Fan J, Yu Y, Lu F, Huang Y.Practical dynamic security regions of bulk power systems.Autom Electr Power Syst 2001; 25(16):6-10.

[17]	Liu Y, Shi X, Xu Y.A hybrid data-driven method for fast approximation of practical dynamic security region boundary of power systems.Int J Electr Power Energy Syst 2020; 117:105658.

[18]	Li B, Wu J, Qiang Z, Qin L, Hao L.Enhanced adaptive assessment on transient stability of power system based on improved deep convolutional generative adversarial network.Autom Electr Power Syst 2022; 46(2):73-82.

[19]	Tan B, Yang J, Lai Q, Xie P, Li J, Xu J.Data augment method for power system transient stability assessment based on improved conditional generative adversarial network.Autom Electr Power Syst 2019; 43(1):149-157.

[20]	Zhou X, Guan X, Sun D, Jiang H, Peng J, Jin Y, et al.Transient stability assessment based on gated graph neural network with imbalanced data in internet of energy.IEEE Internet Things J 2022; 9(12):9320-9331.

[21]	Yang D, Ji M, Zhou B, Bu S, Hu B.Transient stability assessment of power system based on DGL-GAN.Power Syst Technol 2021; 45(8):2934-2945.

[22]	Shi F, Wu J, Wu H, Li B, Ji J, Wang C, et al.Integrated evaluation of power system transient power angle and transient voltage stability margin based on deep learning.Power Grid Technol 2023; 47(2):731-740.

[23]	Tian Y, Wang K, Oluic M, Ghandhari M, Xu J, Li G.Construction of multi-state transient stability boundary based on broad learning.IEEE Trans Power Syst 2021; 36(4):2906-2917.

[24]	Zhan X, Han S, Rong N, Cao Y.A hybrid transfer learning method for transient stability prediction considering sample imbalance.Appl Energy 2023; 333:120573.

[25]	Narasimhamurthi N, Musavi MR.A generalized energy function for transient stability analysis of power system.IEEE Trans Circ Syst 1984; 31(7):637-645.

[26]	Chiang HD, Hirseh MW, Wu FF.Stability regions of nonlinear autonomous dynamical systems.IEEE Trans Automat Contr 1988; 33(1):16-27.

[27]	Yu Y, Liu Y, Qin C, Yang T.Theory and method of power system integrated security region irrelevant to operation states: an introduction.Engineering 2020; 6(7):754-777.

[28]	Feng F, Yu Y.Differential topological characterizations of the dynamic security region of power systems.Proc CSU EPSA 1991; 1(4):48-59.

[29]	Shi X.A hybrid data-driven method for fast approximation of practical dynamic security region boundary of power system [dissertation]. Tianjin University, Tianjin (2020)

[30]	Yu Y, Luan W.A study on dynamic security regions of power systems.Proc CSU EPSA 1990; 2(1):11-21.

RIGHTS & PERMISSIONS

THE AUTHOR

PDF (4940KB)

3481

Accesses

Citation

Detail

Sections

Recommended

Browse

Authors & reviewers

About the journal

Abstract

Graphical abstract

Keywords

Cite this article

1. Introduction

2. Dynamic security region

2.1. Dynamic security region

2.2. Typical differential topological characterizations of DSR

2.3. Practical dynamic security region

3. Space division method

4. Proposed methodology

4.1. Fast search for a small number of critical points

4.2. Critical point augment based on WGAN-GP

4.3. Fit the PDSR boundary

5. Case study

5.1. Setting of the case study

5.2. Case analysis

6. Conclusions

CRediT authorship contribution statement

Declaration of competing interest

Acknowledgments

References

RIGHTS & PERMISSIONS