Identification and Analysis of Cascading Failures in Power Grids with Protective Actions

TL;DR

This paper presents an optimal control-based method to identify initial disturbances that could cause cascading failures in power grids and incorporates protective actions to prevent outages, validated through simulations.

Contribution

It introduces a convex optimization framework that integrates protective actions into the analysis of cascading failures in power grids, providing a novel approach to mitigation.

Findings

Successfully identifies the most disruptive disturbances using the proposed method.

Demonstrates effectiveness of protective actions like load shedding in simulations.

Validates approach on IEEE RTS 24 Bus System.

Abstract

This paper aims to identify and analyze the initial contingencies or disturbances that could lead to the worst-case cascading failures of power grids. An optimal control approach is proposed to determine the most disruptive disturbances on the branch of power transmission system by regarding the disturbances as the control inputs. Moreover, protective actions such as load shedding and generation dispatch are taken into account in a convex optimization framework to prevent the cascading outages of power grids. In theory, the necessary conditions for identifying the most disruptive disturbances are obtained by solving an integrated system of algebraic equations. Finally, numerical simulations are carried out to validate the proposed approach on the IEEE RTS 24 Bus System.

Equations

$$g(P_{ij},c_{ij})=\left\{\begin{array}[]{ll}0,&\hbox{$|P_{ij}|\geq\sqrt{c_{ij}^{2}+\frac{\pi}{2\sigma}}$;}\\ 1,&\hbox{$|P_{ij}|\leq\sqrt{c_{ij}^{2}-\frac{\pi}{2\sigma}}$;}\\ \frac{1-\sin\sigma(P_{ij}^{2}-c_{ij}^{2})}{2},&\hbox{otherwise.}\end{array}\right.$$

$$Y_p^{k+1}=G(P_{ij}^k)\cdot Y_p^k+E_{i_k}{u}_k$$

$$G(P_{ij}^{k})=\left(\begin{array}[]{cccc}g(P_{i_{1}j_{1}},c_{i_{1}j_{1}})&0&.&0\\ 0&g(P_{i_{2}j_{2}},c_{i_{2}j_{2}})&.&0\\ .&.&.&.\\ 0&0&.&g(P_{i_{n}j_{n}},c_{i_{n}j_{n}})\\ \end{array}\right)$$

$$E_{i_k}=diag\left(e_{i_k}^T\right)=diag(\underset{i_k}{\underbrace{0,..,0,1}},0,...,0)\in R^{n\times n}$$

$$P_i=\sum _{j=1}^m{B}_{ij}{\theta }_{ij}=\sum _{j=1}^m{B}_{ij}({\theta }_i-{\theta }_j)$$

$$P=B\theta$$

$$B=\left(\begin{array}[]{cccc}\sum_{i=2}^{m}B_{1i}&-B_{12}&.&B_{1m}\\ -B_{21}&\sum_{i=1,i\neq 2}^{m}B_{2i}&.&B_{2m}\\ .&.&.&.\\ -B_{m1}&-B_{m2}&.&\sum_{i=1}^{m-1}B_{mi}\\ \end{array}\right)$$

$$Y_b^k=A^{T}diag(Y_p^k)A$$

$$P_{ij}^k=e_i^T{Y}_b^k{e}_j(e_i-e_j{)}^T(Y_b^k{)}^{-1^{\ast }}{P}^k,i,j\in I_m=\{1,2,...,m\}$$

$$\begin{array}{rl}& \underset{P^l}{\min }\Vert P^l-P^0{\Vert }^2\\ & s.t.{\underline{P}}_i\le P_i^l\le {\bar{P}}_i\\ & -c_{ij}\le P_{ij}^l\le c_{ij}\end{array}$$

$$\underset{u_k}{\min }J(Y_p^h,u_k)$$

$$J(Y_p^h,u_k)=\Gamma (Y_p^h)+ϵ\sum _{k=0}^{h-1}\frac{\Vert u_k{\Vert }^2}{\max \{0,1-k\}}$$

$$\begin{array}{rl}& 2(P^l-P^0)+\bar{\mu }-\underline{\mu }+\sum _{(i,j)\in \Omega }({\bar{\lambda }}_{ij}-{\underline{\lambda }}_{ij})e_i^T{Y}_b^l{e}_j(Y_b^l{)}^{-1^{\ast }}{e}_{ij}=0\\ & P_i^l-{\bar{P}}_i+{\bar{x}}_i^2=0,{\bar{\mu }}_i(P_i^l-{\bar{P}}_i)=0,{\bar{\mu }}_i-{\bar{z}}_i^2=0\\ & {\underline{P}}_i-P_i^l+{\underline{x}}_i^2=0,{\underline{\mu }}_i(P_i^l-{\underline{P}}_i)=0,{\underline{\mu }}_i-{\underline{z}}_i^2=0\\ & P_{ij}^l-c_{ij}+{\bar{y}}_{ij}^2=0,(P_{ij}^l-c_{ij}){\bar{\lambda }}_{ij}=0,{\bar{\lambda }}_{ij}-{\bar{w}}_{ij}^2=0\\ & P_{ij}^l+c_{ij}-{\underline{y}}_{ij}^2=0,(P_{ij}^l+c_{ij}){\underline{\lambda }}_{ij}=0,{\underline{\lambda }}_{ij}-{\underline{w}}_{ij}^2=0\end{array}$$

$$\begin{array}{rl}& \underset{P^l}{\min }\Vert P^l-P^0{\Vert }^2\\ & s.t.P_i^l-{\bar{P}}_i+{\bar{x}}_i^2=0\\ & {\underline{P}}_i-P_i^l+{\underline{x}}_i^2=0\\ & P_{ij}^l-c_{ij}+{\bar{y}}_{ij}^2=0\\ & P_{ij}^l+c_{ij}-{\underline{y}}_{ij}^2=0\end{array}$$

$$\nabla \Vert P^l-P^0{\Vert }^2+\bar{\mu }-\underline{\mu }+\sum _{(i,j)\in \Omega }({\bar{\lambda }}_{ij}-{\underline{\lambda }}_{ij})e_i^T{Y}_b^l{e}_j(Y_b^l{)}^{-1^{\ast }}{e}_{ij}=0$$

$$\begin{array}{rl}& P_i^l-{\bar{P}}_i+{\bar{x}}_i^2=0,P_{ij}^l-c_{ij}+{\bar{y}}_{ij}^2=0\\ & {\underline{P}}_i-P_i^l+{\underline{x}}_i^2=0,P_{ij}^l+c_{ij}-{\underline{y}}_{ij}^2=0\end{array}$$

$$\begin{array}{rl}& {\bar{\mu }}_i(P_i^l-{\bar{P}}_i)=0,(P_{ij}^l-c_{ij}){\bar{\lambda }}_{ij}=0\\ & {\underline{\mu }}_i(P_i^l-{\underline{P}}_i)=0,(P_{ij}^l+c_{ij}){\underline{\lambda }}_{ij}=0\end{array}$$

$$\begin{array}{rl}& {\bar{\mu }}_i-{\bar{z}}_i^2=0,{\bar{\lambda }}_{ij}-{\bar{w}}_{ij}^2=0\\ & {\underline{\mu }}_i-{\underline{z}}_i^2=0{\underline{\lambda }}_{ij}-{\underline{w}}_{ij}^2=0\end{array}$$

$$\left\{\begin{array}[]{ll}Y_{p}^{k+1}-G(P^{k}_{ij})Y_{p}^{k}-\frac{\max\{0,1-k\}}{2\epsilon}E_{i_{k}}\prod_{s=0}^{h-k-2}\frac{\partial Y_{p}^{h-s}}{\partial Y_{p}^{h-s-1}}\mathbf{1}_{n}=\mathbf{0}\\ \mathcal{F}(Y_{p}^{l},P^{l},\underline{w}_{ij},\bar{w}_{ij},\underline{\lambda}_{ij},\bar{\lambda}_{ij},\underline{y}_{ij},\bar{y}_{ij},\underline{x}_{i},\bar{x}_{i},\underline{z}_{i},\bar{z}_{i},\underline{\mu}_{i},\bar{\mu}_{i})=\mathbf{0}\end{array}\right.$$

$$u_k=\frac{\max \{0,1-k\}}{2ϵ}{E}_{i_k}\prod _{s=0}^{h-k-2}\frac{\partial Y_p^{h-s}}{\partial Y_p^{h-s-1}}{1}_n,k\in I_{h-1}$$

$$Y_p^{k+1}=G(P_{ij}^k)\cdot Y_p^k+E_{i_k}{u}_k$$

$${\left(\frac{\partial Y_p^{k+1}}{\partial u_k}\right)}^T{\lambda }_{k+1}+\frac{ϵ}{\max \{0,1-k\}}\cdot \frac{\partial \Vert u_k{\Vert }^2}{\partial u_k}=0$$

$${\lambda }_k={\left(\frac{\partial Y_p^{k+1}}{\partial Y_p^k}\right)}^T{\lambda }_{k+1}+\frac{ϵ}{\max \{0,1-k\}}\cdot \frac{\partial \Vert u_k{\Vert }^2}{\partial Y_p^k}$$

$$\frac{\partial \mathrm{T}(Y_p^h)}{\partial Y_p^h}-{\lambda }_h=0$$

$$u_k=\frac{\max \{0,1-k\}}{2ϵ}{E}_{i_k}\prod _{s=0}^{h-k-2}\frac{\partial Y_p^{h-s}}{\partial Y_p^{h-s-1}}{1}_n,k\in I_{h-1}$$

$$Y_p^{k+1}-G(P_{ij}^k)Y_p^k-\frac{\max \{0,1-k\}}{2ϵ}{E}_{i_k}\prod _{s=0}^{h-k-2}\frac{\partial Y_p^{h-s}}{\partial Y_p^{h-s-1}}{1}_n=0$$

$$\mathcal{F}(Y_p^l,P^l,{\underline{w}}_{ij},{\bar{w}}_{ij},{\underline{\lambda }}_{ij},{\bar{\lambda }}_{ij},{\underline{y}}_{ij},{\bar{y}}_{ij},{\underline{x}}_i,{\bar{x}}_i,{\underline{z}}_i,{\bar{z}}_i,{\underline{\mu }}_i,{\bar{\mu }}_i)=0$$

$$\gamma =\frac{J(Y_p^h,u)}{J(Y_p^h,0)}.$$

Full text

I Introduction

In the past decades, the world has suffered from several major blackouts such as US-Canada Blackout in 2003 [1], European Blackout in 2006 [2], India Blackout in 2012 and Brazil Blackout in 1999 [3]. All the major blackouts have caused huge economic losses and affected millions of people. Due to the complexity of electrical power systems, it has been a great challenge to understand, analyze and identify the cascading blackouts in practical power grids.

According to the analysis of technological reports, power system blackouts normally go through five stages: precondition, initiating event, cascade events, final state and restoration [3]. The precondition usually happens in the winter or summer peak time due to the excessive power demand. The initiating event (e.g., short-circuit, overload, protection hidden failure, etc.) triggers the chain reaction of branch outages, which starts the cascade events. During the cascade event, power systems may take protective actions such as load shedding and generation dispatch to prevent the cascading failure. If protective actions fail to prevent the further cascades, power systems may end up with the final state of cascading failures, which tends to result in the blackout. As a result, the recovery strategy has to be taken in order to restore the normal state of power grids. To avoid the occurrence of cascading blackouts, it would be desirable to identify the initial malicious disturbances or contingencies before cascading failures, so that the precautions can be taken in advance to eliminate the risk of power system blackouts.

Some identification approaches have been developed to search for the critical branches or initial malicious disturbances that can cause the large-scale disruptions [4, 5, 6, 7, 8, 9, 10, 11]. For instance, some methods are proposed to identify the collections of $n-k$ contingencies based on the event trees [4], line outage distribution factor [5] and other optimization techniques [6, 7, 8]. Nevertheless, these optimization approaches are not efficient to identify the large collections of $n-k$ contingencies that result in cascading blackouts. To address this problem, a “random chemistry” algorithm is designed with the relatively low computational complexity [9]. In addition, an optimal control approach is proposed to identify the initial contingencies by treating these contingencies as the control inputs [10, 11]. The above optimal control approach is able to identify the continuous changes of branch impedance other than direct branch outages as the initial contingencies. Moreover, it can be extended to identify the dangerous fluctuation of injected power on buses caused by the generation of renewable energies or load variations [12]. It is demonstrated that the optimal control approach can effectively determine the worst-case cascades of power grids without protective actions.

As is well known, protective actions such as load shedding and generation dispatch play a key role in preventing cascading blackouts of power grids. This work aims to extend the optimal control approach to identify the malicious contingencies by taking into account protective actions during the cascades. The main challenge in theory is to incorporate protective actions in the framework of optimal control theory, since protective actions are taken according to solutions of a different optimization problem (linear or nonlinear programming). By converting both the optimal control problem and the optimization problem for protective actions into an integrated system of algebraic equations, we manage to formulate and solve the problem of identifying initial contingencies that could lead to the worst-case cascading failures in power grids endowed with protective actions.

The remainder of this paper is organized as follows. Section II formulates the contingency identification problem of power grids involving protective actions in the framework of optimal control theory. Section III presents the theoretical analysis of the optimal control problem. Section IV provides simulation results to validate the proposed approach. Finally, the conclusion is drawn in Section V.

II Problem Formulation

This section presents the problem formulation of identifying the most disruptive disturbances on transmission lines of power grids with $n$ branches and $m$ buses. Figure 1 presents an illustration of identifying the worst-case cascades of power grids with protective actions in the framework of optimal control theory. When the initial contingency or disturbance is added on the branch of power grids, the branch impendence is changed. This leads to the redistribution of power flow on the branches. If the power flow on the branch exceeds its threshold, the branch is severed with the change of network topology. The above branch outage could result in the overloads of other branches and give rise to the further cascades. At the given cascading step (e.g., the $l$-th cascading step in Fig. 1), protective actions (e.g., load shedding and generation dispatch) are taken to prevent the cascading failure. The above cascade process can be described by the state equation in optimal control theory. By designing a cost function to quantify the final disruption level of cascading failures, an optimal control algorithm can be developed to obtain the optimal control inputs, which are exactly the initial contingencies or disturbances that can cause the worst-case cascading failures.

By treating the branch impedance or admittance as the state variable, a state equation can be established to characterize the line outage sequence during the power system cascades. The power flow on each branch is obtained by solving the power flow equation. The protective actions during cascades are implemented at the given cascading step by adjusting the injected power on buses. Finally, an optimal control problem is formulated to search for the initiating event of cascading blackouts with the consideration of load shedding and generation dispatch.

II-A State equation

The state equation is used to describe the line outage sequence of power grids during cascades. Suppose the status of branch or transmission line depends on both the power flow on the branch and its threshold. Thus, an approximate function on the branch status is designed to characterize the operation of circuit breaker in protective relays as follows

[TABLE]

where $\sigma$ is a tunable positive parameter. $P_{ij}$ refers to the power flow on the branch between Bus $i$ and Bus $j$, and $c_{ij}$ denotes its threshold of power flow. $P_{ij}$ can be obtained by solving the power flow equation. Notably, the function $g$ is differentiable with respect to $P_{ij}$, and it approximates to the step function as the parameter $\sigma$ increases. It can reflect the physical characteristic of protective relays while contributing to theoretical analysis on the identification of initial contingencies using optimal control theory.

The line outage sequence can be described by the following state equation (i.e., branch outage model in Figure 1)

[TABLE]

where $Y_{p}^{k}=(y^{k}_{p,1},y^{k}_{p,2},...,y^{k}_{p,n})^{T}\in R^{n}$ refers to the vector of branch admittance, and the control input or disturbance is denoted by $u_{k}=(u_{k,1},u_{k,2},...,u_{k,n})^{T}$. In addition, $G(P_{ij}^{k})$ is a diagonal matrix

[TABLE]

and it describes the status of each branch at the $k$-th cascading step. The matrix $E_{i_{k}}$ is constructed as

[TABLE]

and it enables us to select the $i_{k}$-th branch for adding the initial disturbances.

II-B DC power flow equation

For high-voltage transmission systems, the DC power flow is a good substitute for the AC power flow [13], and it is computationally efficient and immune to the numerical non-convergence. Thus, the DC power flow equation is employed in this work. The classic DC power flow equation is given by

[TABLE]

where $P_{i}$ and $\theta_{i}$ refer to the injection power and voltage phase angle of Bus $i$, respectively. $B_{ij}$ represents the mutual susceptance between Bus $i$ and Bus $j$. The symbol $m$ denotes the number of buses in power networks. Equation (2) can be rewritten in matrix form [14]

[TABLE]

where $P=(P_{1},P_{2},...,P_{m})$, $\theta=(\theta_{1},\theta_{2},...,\theta_{m})$ and

[TABLE]

It is worth pointing out that $B$ can be obtained by removing the real part of $Y_{b}$. In fact, the nodal admittance matrix $Y^{k}_{b}$ at the $k$-th cascading step can be obtained as

[TABLE]

where $A$ denotes the branch-bus incidence matrix [15]. Thus, the power flow on the branch between Bus $i$ and Bus $j$ can be computed as

[TABLE]

where $e_{i}$ is the $m$-dimensional unit vector with $1$ at the $i$-th position and [math] elsewhere. $P^{k}$ refers to the vector of injected power on buses at the $k$-th cascading step. The symbol $-1^{*}$ represents a generalized inverse for solving DC power flow equation [10].

Remark II.1**.**

During the cascading blackout, power network may be divided into several subnetworks ($i.e.$, islands), which can be identified by analyzing the nodal admittance matrix $Y^{k}_{b}$. To solve the DC power flow equation, the generator bus connected to the largest generating station is selected as the new slack bus in the subnetwork. And thus the power variation of slack bus accounts for a small percentage of its generating capacity. If there is no generator bus in the subnetwork, the power flow is zero on each branch of this subnetworks.

II-C Protective actions

If generation dispatch and load shedding are taken into account in the formulation, $P^{k}$ has to be updated at certain steps of cascading failure. For simplicity, suppose that load shedding and generation control are implemented at the $l$-th cascading step ($1<l<h$). This implies that $P^{k}=P^{l}$ for $k\geq l$. Thus, a nonlinear programming problem can be proposed to allow for load shedding and generation dispatch as follows

[TABLE]

where $P^{l}=(P_{1}^{l},P_{2}^{l},...,P_{m}^{l})^{T}$, and $P^{0}$ denotes the vector of original injected power on buses. The symbols $\underline{P}_{i}$ and $\bar{P}_{i}$ denote the upper and lower bounds of injected power on Bus $i$, respectively. The cost function in (3) quantifies the changes of injected power on buses due to load shedding and generation control. Essentially, the objective of Optimization Problem (3) is to achieve the minimum adjustment of injected power on buses while preventing the further branch outages of power grids.

Remark II.2**.**

The linear programming formulation can also be adopted to allow for protection actions in power systems [16]. The proposed approach is also applied to the linear programming formulation. In addition, it also has a chance to be extended for dealing with protection actions at multiple cascading steps.

II-D Cost function

Next, a cost function of optimal control problem is presented to quantify the disruption level of power grids at the end of cascading failures. Suppose the cascades come to an end at the $h$-th cascading step, and the cost function is designed as

[TABLE]

with

[TABLE]

and the adjustment of injected power on buses $P^{k}$ for the system protection is implemented according to the solutions to Optimization Problem (3) at the $l$-th cascading step. The first term in the cost function (i.e., $\Gamma(Y_{p}^{h})$) describes the final status of cascading failures (e.g., network connectivity or power flow), and the second term characterizes the accumulated control energy. In addition, $\epsilon$ is a positive weight, and the symbol $\|\cdot\|$ represents the $2$-norm. Normally, $\epsilon$ is small enough so that more efforts are taken to minimize the first term in the cost function.

Remark II.3**.**

The second term in the cost function ensures that only initial disturbances can be added on the branch of power transmission system. This is due to $\max\{0,1-k\}=1$ for $k=0$ and $\max\{0,1-k\}=0$ for $k\geq 1$. When $\max\{0,1-k\}=0$, the second term goes to the infinity. This implies that it is impossible to minimize the cost function when $k\geq 1$.

III Theoretical Results

This section provides theoretical results on the identification of initial contingencies on branches that can cause the catastrophic cascading failures of power grids. First of all, an equivalent condition is presented for Problem (3).

Lemma III.1**.**

The optimal solutions of Optimization Problem (3) are equivalent to the solutions of the following system of algebraic equations

[TABLE]

where $Y^{l}_{b}=A^{T}diag(Y^{l}_{p})A$, $e_{ij}=e_{i}-e_{j}$, $i\in I_{m}$, $(i,j)\in\Omega$, $\bar{\mu}=(\bar{\mu}_{1},\bar{\mu}_{2},...,\bar{\mu}_{m})^{T}$ and $\underline{\mu}=(\underline{\mu}_{1},\underline{\mu}_{2},...,\underline{\mu}_{m})^{T}$. And the symbol $\Omega$ denotes the set of branches in power systems with the cardinality $|\Omega|=n$ (i.e., the number of branches).

Proof.

The inequality constraints in the optimization problem (3) can be converted into the equality constraints by introducing the unknown variables $\bar{x}_{i}$, $\underline{x}_{i}$, $\bar{y}_{ij}$ and $\underline{y}_{ij}$ as follows.

[TABLE]

According to the Karush-Kuhn-Tucker (KKT) conditions, the necessary condition for the solutions to Optimization Problem (6) can be obtained as follows. Specifically, we have

[TABLE]

with $\nabla\|P^{l}-P^{0}\|^{2}=2(P^{l}-P^{0})$ for the stationarity condition and the constraints

[TABLE]

The conditions for complementary slackness are given by

[TABLE]

Moreover, the dual feasibility can be described by

[TABLE]

where $e_{ij}=e_{i}-e_{j}$, $i\in I_{m}$, $(i,j)\in\Omega$, $\bar{\mu}=(\bar{\mu}_{1},\bar{\mu}_{2},...,\bar{\mu}_{m})^{T}$ and $\underline{\mu}=(\underline{\mu}_{1},\underline{\mu}_{2},...,\underline{\mu}_{m})^{T}$. Since the cost function in Optimization Problem (3) is a convex function and the inequality constraints are affine, the above necessary conditions are also sufficient for optimality. This implies that the solutions to Optimization Problem (3) are equivalent to solutions to the system of algebraic equations (6). The proof of this lemma is thus completed. ∎

Remark III.1**.**

As we can see, System (5) is composed of $(7m+6n)$ equations and $(7m+6n)$ additional unknown variables (i.e., $P^{l}$, $\underline{w}_{ij}$, $\bar{w}_{ij}$, $\underline{\lambda}_{ij}$, $\bar{\lambda}_{ij}$, $\underline{y}_{ij}$, $\bar{y}_{ij}$, $\underline{x}_{i}$, $\bar{x}_{i}$, $\underline{z}_{i}$, $\bar{z}_{i}$, $\underline{\mu}_{i}$, $\bar{\mu}_{i}$). Note that $Y^{l}_{b}$ contains the existing unknown variables $Y_{p}^{l}$ in the state equation (1).

The equivalent conditions in Lemma III.1 allows to obtain the necessary conditions for the optimal control problem (4) as follows

Theorem III.1**.**

The necessary conditions for the optimal control problem (4) with protective actions according to (3) are given by solving the system of algebraic equations as follows

[TABLE]

where the second equation represents System (5), and the optimal adjustment of injected power on buses for protective actions satisfies $P^{k}=P^{l}$ for $k\geq l$, and $P^{l}$ is the solution to Problem (3). In addition, the optimal control input is given by

[TABLE]

Proof.

With Pontryagin’s maximum principle in optimal control theory for the discrete-time system [17], the necessary conditions for the optimal control problem (4) can be determined as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\mathbf{0}=(0,0,...,0)^{T}\in R^{n}$. By reorganizing the above equations, we can obtain the optimal control input

[TABLE]

and the system of algebraic equations

[TABLE]

where the vector $P^{l}$ is determined by the solutions to System (5), which can be rewritten as

[TABLE]

for ease of notation according to Lemma III.1. By combining Equations (9) and (10), an extended system of algebraic equations is obtained with $7m+(6+h)n$ equations and $7m+(6+h)n$ unknown variables. By substituting the solutions of this extended system into (8), we can identify the most disruptive disturbances for the cascades of power grids with protective actions at the given cascading step. This completes the proof ∎

Remark III.2**.**

The extended system of algebraic equations (7) can be solved using the numerical solver in numerical-analysis softwares. In this way, load shedding and generation dispatch can be taken into account in the proposed optimal control formulation of identifying the worst-case cascading failures.

Remark III.3**.**

Since the solutions to System (7) can only provide necessary conditions for the optimal control problem (4), extensive numerical simulations have to be conducted to search for the feasible solutions in practice. By comparing the values of cost function with the identified initial disturbances (i.e., control inputs), it is expected to obtain the optimal or suboptimal solutions to Problem (4) if numerical simulations are repeated for a sufficiently large number of times.

IV Simulation and Validation

To demonstrate the effectiveness of the proposed identification approach based on optimal control, numerical simulations are conducted on the IEEE RTS 24 Bus System to determine the initial disruptive disturbances on each branch (see Fig. 2). In addition, we make a comparison on simulation results between the cases with and without protective actions (i.e., load shedding and generation dispatch), respectively.

IV-A Parameter setting

Per unit values are adopted with the base value of power $100$ MVA. Moreover, the solver “fsolve” in Matlab is employed to solve the integrated system of algebraic equations (7). Other parameters are given as follows: $\sigma=5\times 10^{4}$, $\epsilon=10^{-4}$, $h=10$, and $\Gamma(Y_{p}^{h})=\|Y_{p}^{h}\|^{2}/2$ in the cost function (4). For each branch, numerical simulations are carried out to solve System (7) for $10$ times, and the worst-case disturbances are selected by comparing the values of cost function. A performance index $\gamma$ is defined to quantify the disruption level of cascades as follows

[TABLE]

Intuitively, the index $\gamma$ is the ratio between the final cost of power grids with the control input $u$ and that without any control inputs, and a smaller $\gamma$ indicates a worse cascade of power systems. The power flow threshold on each branch is $10\%$ larger than the normal power flow on the corresponding branch of power systems without any disturbances.

IV-B Validation and comparison

In the simulations, generation dispatch and load shedding are implemented at the 4-th cascading step (i.e., $l=4$) according to solutions of Optimization Problem (3). Figure 3 presents the initial disturbances (i.e., control inputs) on each branch identified by the proposed optimal control approach and the resulting normalized costs (i.e., the index $\gamma$). Specifically, the blue bars denote the control inputs and normalized costs without generation dispatch and load shedding, while the green bars represent those with generation dispatch and load shedding. As we can observe in the upper panel of Figure 3, the height of green bar is not smaller than that of blue bar for each branch. This indicates that the larger initial disturbances (i.e., the magnitude of control input) are required to trigger the worst-case cascades of power grids with generation dispatch and load shedding compared to those without generation dispatch and load shedding. This demonstrates that protection actions are able to effectively enhance the robustness of power systems and relieve the final disruption level after malicious disturbances. It is worth pointing out that the branches with equal heights of blue bar and green bar (e.g., Branch 1, Branch 2, Branch 3, Branch 4, Branch 5, Branch 6, etc) are directly severed by the initial disturbances (i.e., control inputs). The lower panel of Figure 3 demonstrates that the cascades with generation dispatch and load shedding are less disruptive on the whole except for the cascades triggered by the initial disturbances on Branch 28, Branch 32 and Branch 33. This is because the optimal solutions for generation dispatch and load shedding (i.e., solutions to the optimization problem (3)) are not obtained by the numeric solver at the specified cascading step. Thus, the adjustment of injected power on buses (i.e., generation dispatch and load shedding) actually deteriorates branch overloads and results in the worse disruptions of cascades in the end.

V Conclusions

In this paper, we investigated the problem of identifying the initial contingencies that result in the worst-case cascading failures of power grids. In particular, power grids are equipped with protection devices to prevent the cascading blackouts by load shedding or generation dispatch in time. Moreover, a theoretical framework was proposed to allow for both the identification of the most disruptive disturbances and the optimal adjustment of injected power buses for the protection of power girds. Numerical simulations were conducted to better understand the effect of protection actions on the identification of initial malicious contingencies and the final disruptions of cascading failures. In this work, the deterministic cascading failure paths are taken into account for the identification of initial disruptive disturbances and the implementation of protective actions. In practice, power grids are subject to uncertainties (e.g., hidden failure, device aging, human errors, etc). Therefore, we will consider the contingency identification of power grids with uncertain cascading failure paths in the next step.

Acknowledgment

This work is partially supported by the Future Resilience System Project at the Singapore-ETH Centre (SEC), which is funded by the National Research Foundation of Singapore (NRF) under its Campus for Research Excellence and Technological Enterprise (CREATE) program. It is also supported by Ministry of Education of Singapore under Contract MOE2016-T2-1-119.