A Two-Stage Optimal Preventive Control Model Incorporating Transient Stability Constraints in the Presence of Multi-Resource Uncertainties

Abstract

The volatility and uncertainty introduced by increasingly integrated renewable energy pose challenges to the reliable and stable operation of the power system. To mitigate the operation risks, a two-stage optimal preventive control model that incorporates transient stability constraints and considers uncertainties from multiple resources is proposed. First, the uncertainties of different re-sources are modeled, with which the non-sequential Monte Carlo sampling method is used to correspondingly generate the scenarios. Thereafter, a two-stage control model that balances operational safety and economy and realizes preventive control and emergency control is built. The operation schedule from the preventive control stage aims to minimize the transient stability probability and operation costs. If any faults destabilize the system, the emergency control stage will be activated immediately to help the system recover stability with minimal control costs. To expedite the solving of the two-stage model, a multi-objective particle swarm algorithm based on entropy-TOPSIS is proposed. Finally, the effectiveness of the proposed model and solving algorithm are validated with the modified IEEE118 node system.

1. Introduction

With the increasing integration of renewable energy that introduces substantial uncertainty, the stable operation of the system is facing severe challenges [1,2]. Traditional optimal operation models usually ignore system transient stability, resulting in the generated operation schedules failing to satisfy the system’s stable operation requirements. Therefore, it is necessary to consider the transient stability constraint in the optimal operation model and thereby form an optimal preventive control model. The transient stability-constrained optimal operation problem is a topic of much recent interest as it offers a compromise between economic efficiency and safe operation of power systems. The transient stability-constrained model uses differential-algebraic equations (DAEs) to describe system dynamics and transient stability constraints [3] and is therefore usually cast as a semi-infinite programming-based model, which re-quires a significant computational effort to solve. Therefore, it is necessary to propose a new optimal control method which can take the transient stability into consideration and can be easily solved.

Considering the involved uncertainty, paving the way to the optimal preventive control of the systems should resolve three problems: (1) modeling the uncertainty; (2) building the preventive control model; and (3) solving the optimal preventive control model.

• Modeling the uncertainty: Reference [4] focuses on modeling the uncertainty of wind and solar energy, mainly by adopting the Burr distribution model [5], Weibull distribution model [6], Gamma distribution model [7], etc. The Weibull distribution model has been widely used due to its simple principles and convenience in calculation. In addition, the structure of a power grid can also be impacted by unexpected faults. The uncertainty of the occurrence of faults is analyzed from three aspects: fault type, fault location, and recovery time. In existing works, the uncertainty of fault types is described with discrete probability distribution models [8], and the uncertainty of fault recovery time is usually described with a normal distribution [9]. Similarly, the load uncertainty, mainly caused by the fluctuation of active power injected into the node, can be described by a normal distribution [10]. In the literature, random samples are usually generated via Monte Carlo sampling simulations [11] and Latin hypercube sampling simulations [12], and the generated scenarios are reduced using clustering methods such as the k-means method.
• Building preventive control model: Non-anticipative changes on grid components and the aroused redistribution of currents can destabilize the system. Preventive control prevents possible system instability by adjusting the current operation status with the cost of a higher operating cost. References [13,14,15,16,17] studied preventive control strategies on hourly time basis. The difference between these models is mainly their objectives, which can be roughly divided into three categories: eco-nomic cost minimization [13]; reliability maximization [14,15]; and joint economy and reliability [16,17]. Reference [18] mainly studies the online preventive control and proposes a method for safety and stability control considering large-scale wind energy. However, preventive control is indeed a multi-timescale problem that is solved at different time scales step by step. Therefore, it is reasonable to build the prevention and control model in two parts, namely, preventive control and emergency control.
• Solving methodology: A optimization model can be cast as a given programming problem, such as mixed integer linear programming (MILP), and then solved with commercial solvers. On one hand, this methodology usually requires linearization of the original model which is originally nonlinear and nonconvex; on the other hand, for large-scale problems, the computational burden is heavy and sometimes unacceptable. Intelligent algorithms, such as particle swarm optimization algorithms [19] and grey wolf algorithms [20], are wildly used in solving complicated optimization model in recent years. The advantage is that they can be applied to different models without additional processing, but their solution stability is expected to be further improved.

This paper proposes a two-stage optimal preventive control model that considers transient stability constraints and uncertainties from multiple resources. First, the un-certainty from the generation units, the grid, and the loads are modeled, and a set of scenarios is generated via the Monte Carlo sampling method. Second, a two-stage control model that balances system operation safety and economy is proposed. The first stage is the preventive control stage, with which the transient stability probability and operating costs are considered based on the faults in uncertain scenarios. The second stage is the emergency control strategy, which is executed in the event of a fault to minimize the control cost and recover system stability. To improve efficiency, a multi objective particle swarm algorithm based on entropy-TOPSIS is proposed to solve the problem. Finally, the effectiveness of the proposed methodology is demonstrated with the modified IEEE118 node system. The main contributions of the manuscript are as follows:

(1).

The paper proposes a two-stage optimal preventive control model that incorporates transient stability constraints and addresses the challenges aroused by the volatility of renewable energy resources.

(2).

The non-sequential Monte Carlo-based scenario generation method is used to simulate the uncertainties introduced by multiple resources. This provides a complete under-standing of the system’s behavior and enhances the reliability of decisions in grid operation.

(3).

An improved multi-objective particle swarm algorithm (IMOSPO) based on entropy-TOPSIS is proposed to efficiently solve the proposed model.

The rest of this manuscript is organized as follows: Section 2 gives the uncertainty models for multiple resources; Section 3 describes the details of the proposed two-stage optimal preventive control model; the solution methodology is introduced in Section 4; the case study is presented in Section 5; finally, Section 6 concludes this paper.

2. Probability-Based Uncertainty Models

2.1. Modeling the Uncertainty of System Components

2.1.1. Modeling the Uncertainty of Generation Units

With the sustainable development of renewable energy, the proportion of renew-able generation with substantial uncertainty on the generation side gradually increases. Renewable resources have become the source of uncertainty on the generation side. Especially for the wind turbine, the uncertainty of wind speed and turbine failure may frequently result in wind power output fluctuation. Thus, in this paper, the wind tur-bine is taken as an example to introduce the uncertainty modeling of the generation side. In addition, the uncertainty of other renewable energy resources, such as PV, can be described following the model of the wind turbine.

Weibull distribution is widely used to describe wind speed distribution. The prob-ability density function of Weibull distribution can be written as in (1):

$$f\left(v\right)=\frac{k}{c}·{\left(\frac{v}{c}\right)}^{k-1}·exp\left[-{\left(\frac{v}{c}\right)}^k\right]$$

(1)

where c is the size factor; k is the shape factor; and v is the input wind speed.

The relationship of the output power $P^W$ and input wind speed of a wind turbine $v_i$ can be expressed as in (2):

$$P^W=\left\{\begin{array}{c}0,v_i>v_o \mathrm{or} v_i<v_{in}\\ P_N^W·\frac{\left(v_i-v_{in}\right)}{\left(v_N-v_{in}\right)}, v_{in}\le v_i\le v_N\\ P_N, v_N\le v_i\le v_o \end{array}\right.$$

(2)

2.1.2. Modeling the Uncertainty of the Grid

The uncertainty faced by the power grid is mainly aroused by the unexpected failures on transmission lines and thus can be featured by the fault type, fault location, and the corresponding recovery time. Therefore, describing these three features is im-portant to model the uncertainty of the power grid.

Probability Distribution of the Fault Type

According to the statistical results of electrical faults that happened in North America from the Institute of Electrical and Electronics Engineers (IEEE), the discrete probabilities of different fault types are shown in

Table 1. Probabilities of different fault types.

Fault Type	Three-Phase	Two-Phase	Single-Phase	Phase-to-Phase	Total
Statistical results %	1	2	93	4	100

.

Probability Distribution of the Fault Location

The fault location is described by the distance between the fault location and the end of the line, and its probability distribution is modeled with a discrete distribution as shown as in (3):

$$f\left(i\right)=\left\{\begin{array}{c}0.2, \frac{L_i}{L}<20\%\\ 0.6, 20\%<\frac{L_i}{L}<80\%\\ 0.2, \frac{L_i}{L}>80\%\end{array}\right.$$

(3)

where $L_i$ is the distance between the fault location and the end of a line; and $L$ represents the total length of the line.

Probability Distribution of the Recovery Time

The process of fault recovery on the grid cannot be ignored because the fault recovery time will reflect the time duration that the impact of the fault lasts on the grid. Considering the different operating statuses of relay protection systems and communication networks, the fault recovery time is not constant and is modeled with normal distribution as shown in (4) [21]:

$$f\left(t_i\right)=\frac{1}{\sqrt{2\pi }{\sigma }_t}exp\left[-\frac{{(t_i-{\mu }_t)}^2}{2{\sigma }_t{}^2}\right]$$

(4)

where the $t_i$ is the fault recovery time; ${\mu }_t$ is the expectation of the fault recovery time; and ${\sigma }_t$ is the variance of the fault recovery time.

2.1.3. Uncertainty Modeling of Load

Load is randomness and volatility, whose distribution can be described by normal distribution:

$$f\left(P_n\right)=\frac{1}{\sqrt{2\pi }{\sigma }_n}exp\left[-\frac{{(P_n-{\mu }_n)}^2}{2{\sigma }_n{}^2}\right]$$

(5)

where $P_n$ is the active power, ${\mu }_n$ is the expectation of active power, and ${\sigma }_n$ is the variance of active power.

2.2. Non-Sequential Monte Carlo Based Scenarios Generation Method

The Monte Carlo method is often used in analyzing the reliability and stability of large complex networks. Based on the law of large numbers, the core idea of the Monte Carlo method is to use enough repeated random experiments to solve problems with probabilistic interpretations and obtain statistical results. Furthermore, the Monte Carlo method can be classified as the category of non-sequential Monte Carlo methods and the category of sequential Monte Carlo methods.

The non-sequential Monte Carlo method uses the random sampling method to generate the operation status of each component. It then analyzes the operation status of the entire system according to the operation statuses of its components. Evidently, the uncertainties contained by the components are finally aggregated as the uncertainty of the entire system. The non-sequential Monte Carlo method is simple and computationally friendly. Thus, it is ideal for analyzing the reliability of large-scale power systems with speed requirements.

To generate a stochastic scenario of the entire system using the non-sequential Monte Carlo method, the operation status of each component needs to be sampled first. Considering that the system contains $N$ components, $e_k$ is used to indicate the component $k$’s operation status, where $k=1,2,3,\dots ,N$. $e_k=1$ represents equipment $k$ operating normally; $e_k=0$represents equipment $k$ operating abnormally. The samples of every component can be written as (6):

$$e_k=\left\{\begin{array}{c}1 T_k\le H_k\\ 0 0\le H_k\le T_k\end{array}\right.$$

(6)

where $T_k$ is the probability of component k being abnormal; and$H_k$ is a random number that obeys uniform distribution in the interval of [0, 1].

The above steps can be repeated to generate the operation statuses of the remaining components. When the statuses of all components are determined, the random scenario generation of the entire system is completed.

Considering the operation status of the simulated system is $E_i$, the corresponding probability of the operation status of the system can be indicated by $P\left(E_i\right)$ as in (7):

$$P\left(E_i\right)={{\displaystyle \prod }}_{k=1}^{N_o}{T}_k{{\displaystyle \prod }}_{k=1}^{N_o-N}{H}_k$$

(7)

where $N_o$ is the number of abnormal components in the system. Based on the law of large numbers, when the number of sampling times goes to infinite, the average value of samples will be infinitely close to the expected value of samples.

3. The Two-Stage Optimal Preventive Control Model

The preventive control model proposed in this paper is a two-stage optimization model. The first stage focuses on prevention control and considers the power and re-serve from generation units and energy storage systems as variables. In this stage, the goal is to minimize the total cost and the probability of system transient instability while considering the economic and stable aspects of the system. The second stage is the emergency response stage. With a fault event, evaluation is made to determine if the system will experience any instability. If instability is detected, the emergency control stage will be triggered. During the second stage, the decision variables involve the usage of reserves from generation units and energy storage systems. The focus of this stage is to minimize the costs of delivering reserves and to achieve a cost-effective restoration of the system’s transient stability.

3.1. Model of Preventive Control Stage

3.1.1. Objective Function

The preventive control stage needs to balance safety and the economy. Thus, the objective of the proposed model can be written as in (8), where $f_1$ represents the total cost; and $f_2$ stands for the stability of the system. The cost of the system consists of the expected operation costs of the generation units. The cost function $f\left(p_{t_1,0}\right)$ and $f\left(p_{t,s}\right)$, respectively, represent the total cost of the system at the time $t_1$ in the base scenario (scenario $0$) and at time $t$ in scenario $s$. The stability of the system can be reflected by the probability of transient stability, which is represented by (13). The stability of the system mainly depends on the fluctuation of power and the failures of the components. $\beta$ is the transient stability probability.

$$\mathrm{Min} \mathrm{F}\left(x\right)=\left\{f_{1,}{f}_2\right\}$$

(8)

$$f_1=f\left(p_{t_1,0}\right)+{{\displaystyle \sum }}_{s=1}^{n_s}{\pi }_s{{\displaystyle \sum }}_{t=t_1+1}^{n_T}f\left(p_{t,s}\right)$$

(9)

$$f\left(p_{t_1,0}\right)={{\displaystyle \sum }}_{i=1}^{n_G}[(a_i{\left(p_{i,t_1,0}^G\right)}^2+b_i{p}_{i,t_1,0}^G+c_i+r_{i,t_1}^{G,UP}{C}^{UP}+r_{i,t_1}^{G,DN}{C}^{DN}\left]+{{\displaystyle \sum }}_{e=1}^{n_s}\left(-p_{e,t_1,0}^C{C}_{e,t1}^C+p_{e,t_1,0}^D{C}_{e,t1}^D\right)\right]$$

(10)

$$f\left(p_{t,s}\right)={{\displaystyle \sum }}_{i=1}^{n_G}[(a_i{\left(p_{i,t,s}^G\right)}^2+b_i{p}_{i,t,s}^G+c_i+r_{i,t}^{G,UP}{C}^{UP}+r_{i,t}^{G,DN}{C}^{DN}\left]+{{\displaystyle \sum }}_{i=1}^{n_s}\left(-p_{e,t,s}^C{C}_{e,t}^C+p_{e,t,s}^D{C}_{e,t}^D\right)\right]$$

(11)

$$p_{t,s}=\left\{p_{i,t,s}^G,p_{e,t,s}^C,p_{e,t,s}^D,r_{i,t}^{G,UP},r_{i,t}^{G,DN}\right\}$$

(12)

$$f_2=\beta =1-\frac{n_{ss}}{S}$$

(13)

In (10)–(13), $p_{t,s}$ represents the vector of decision variables; $a_i$, $b_i$, and $c_i$ are, respectively, the cost coefficients of different sets of generation units; ${\pi }_s$ is the probability of uncertainty scenario $s$; $n_G$ is the total number of generation units; $n_T$ is the number of time intervals of the preventive control stage; $n_G$ is the number of generators; $n_s$ is the number of storage devices; $n_{ss}$ is the number of scenarios that the system is in a statable status; $p_{i,t,s}^G$ is the active power output of the $i^{th}$ generation unit at time $t$ in scenario $s$; $p_{e,t,s}^C$ and $p_{e,t,s}^D$ are the charging and discharging powers, respectively, of energy storage system $e$; $C_{e,t}^C$ and $C_{e,t}^D$ are the charging and discharging costs, respectively, of energy storage system $e$; $r_{i,t}^{G,UP}$ and $r_{i,t}^{G,DN}$ represent the up and down reserve capacity from generation units, respectively; $C^{UP}$ and $C^{DN}$ represent, respectively, the unit costs for the up and down reserve capacity from generation units; the probability of the system’s transient instability is $\beta$; and $S$ represents the number of scenarios.

3.1.2. Constraints

The constraints could be divided into five parts: system operation constraints; constraints of energy storage; constraints of generation unit; transient process constraints; and stability margin constraints.

(1)

System operation constraints

Constraints (14) and (15), respectively, represent the active and reactive power balance of the system level. The probabilistic constraints of node voltage limitation and transmission line capacity are shown in (16) and (17), respectively.

$$p_{b,t,s}^G+p_{b,t,s}^W-p_{b,t,s}^C+p_{b,t,s}^D-p_{b,t,s}^L-u_{b,t,s}{{\displaystyle \sum }}_{j=1}^{n_b}{u}_{j,t,s}\left(G_{bj}cos{\delta }_{bj,t,s}+B_{bj}sin{\delta }_{bj,t,s}\right)=0$$

(14)

$$q_{b,t,s}^G+q_{b,t,s}^W-q_{b,t,s}^L+u_{b,t,s}{{\displaystyle \sum }}_{j=1}^{n_b}{u}_{j,t,s}\left(B_{bj}cos{\delta }_{bj,t,s}-G_{bj}sin{\delta }_{bj,t,s}\right)=0$$

(15)

$$Pro\left\{{\underset{\_}{U}}_b\le u_{b,t,s}\le {\overline{U}}_b\right\}>{\alpha }_U$$

(16)

$$Pro\left\{-{\overline{S}}_l\le s_{l,t,s}\le {\overline{S}}_l\right\}>{\alpha }_S$$

(17)

In (14)–(17), $n_b$ denotes the number of nodes; $p_{b,t,s}^W$ and $q_{b,t,s}^W$ are the active and reactive power output of the wind turbines at node $b$ in the scenario $s$, respectively; $p_{b,t,s}^L$ and $q_{b,t,s}^L$, respectively, denote the active and reactive power of the transmission line; $q_{b,t,s}^G$ is the reactive power of the $i^{th}$ generation unit at time t in scenario $s$; $G_{bj}$ and $B_{bj}$ are the conductivity and admittance of branch $bj$, respectively; $u_{b,t,s}$ and $u_{j,t,s}$ are, respectively, the voltage amplitude of node $b$ and $j$ in scenario $s$; ${\delta }_{bj,t,s}$ represents the difference of phase angle between node $b$ and $j$ in scenario $s$; ${\overline{U}}_b$ and ${\underset{\_}{U}}_b$ are, respectively, the upper and lower limits of the voltage amplitude of node $i$; $s_{l,t,s}$ and ${\overline{S}}_l$ are, respectively, the apparent power of the $i^{th}$ line and the corresponding upper limit; ${\alpha }_U$ and ${\alpha }_S$ are, respectively, the probability thresholds of node voltage and apparent power limitations.

(2)

Constraints of energy storage

Energy storage’s charging and discharging power limits are enforced with (18) and (19). The upward reserve of energy storage can be provided by discharging more and the downward reserve can be provided by charging more. Constraint (20) enforces that both the upward reserve and downward reserve should be positive. The energy evolution function is shown in constraint (21), which builds the relationship between the $e_{e,t-1}^S$ and $e_{e,t}^S$. Constraint (22) enforces the upper and lower bounds of stored energy. Moreover, constraint (23) indicates that simultaneously charging and discharging is not allowed.

$${\underset{\_}{P}}_e^C{c}_{e,t}\le p_{e,t,s}^C+r_{e,t}^C\le {\overline{P}}_e^C{c}_{e,t}$$

(18)

$${\underset{\_}{P}}_e^D{d}_{e,t}\le p_{i,t,s}^D+r_{e,t}^D\le {\overline{P}}_e^D{d}_{e,t}$$

(19)

$$r_{e,t}^C\ge 0;r_{e,t}^D\ge 0$$

(20)

$$e_{e,t}^S=e_{e,t-1}^S+p_{e,t}^C{\eta }_e^C\Delta t-p_{e,t}^D/{\eta }_e^D\Delta t$$

(21)

$${\underset{\_}{E}}_e^S\le e_{e,t}^S\le {\overline{E}}_e^S$$

(22)

$$c_{e,t}+d_{e,t}\le 1$$

(23)

In (18)–(23), ${\overline{P}}_e^C$ and ${\underset{\_}{P}}_e^C$ are, respectively, the upper and lower bounds of charging power; $c_{e,t}$ and $d_{e,t}$ are, respectively, the charging and discharging status of energy storage system $e$; $\Delta t$ is the length of the time interval; $E_e^{\mathrm{S}}$ represents the stored energy; ${\overline{E}}_e^S$ and ${\underset{\_}{E}}_e^S$ are upper and lower bounds of stored energy system $e$, respectively; ${\eta }_e^C$ and ${\eta }_e^D$ are the charging and discharging efficiencies, respectively.

(3)

Constraints of generation units

Several aspects should be considered for the operation of generation units. To provide the upward reserve and downward reserve, the upper and lower bound constraints of generation units are modified via (24). The ramping constraint shown in (25) limits the output change. The total upward reserve and downward reserve from generation units and energy storage systems should satisfy the minimum reserve requirement of the system as in (26) and (27). The reactive power outputs of generation units are limited by (28). The active power outputs of wind turbines should be less than the corresponding forecasts as in (29). Similarly, the reactive power outputs of wind turbines are limited by (30).

$${\underset{\_}{P}}_i^G+r_{i,t}^{G,DN}\le p_{i,t,s}^G\le {\overline{P}}_i^G-r_{i,t}^{G,UP}$$

(24)

$${\underset{¯}{RP}}_i^G\le p_{i,t,s}^G-p_{i,t-1,s}^G\le {\overline{RP}}_i^G$$

(25)

$${{\displaystyle \sum }}_{i=1}^{n_G}{r}_{i,t}^{G,UP}+{{\displaystyle \sum }}_{e=1}^{n_S}{r}_{e,t}^D\ge R_t^{UP}$$

(26)

$${{\displaystyle \sum }}_{i=1}^{n_G}{r}_{i,t}^{G,DN}+{{\displaystyle \sum }}_{e=1}^{n_S}{r}_{e,t}^C\ge R_t^{DN}$$

(27)

$${\underset{\_}{Q}}_i^G\le q_{i,t,s}^G\le {\overline{Q}}_i^G$$

(28)

$$0\le p_{b,t,s}^W\le {\overline{P}}_{b,t,s}^W$$

(29)

$$0\le q_{b,t,s}^W\le {\overline{Q}}_{b,t,s}^W$$

(30)

In (24)–(30), ${\overline{P}}_i^G$ and ${\underset{\_}{P}}_i^G$ are, respectively, the upper and lower bounds of the power outputs of generation units; ${\overline{RP}}_i^G$ and ${\underset{¯}{RP}}_i^G$, respectively, denote the upper and lower bounds of the ramping capacity of generation units; $R_t^{UP}$ and $R_t^{DN}$ are, respectively, the minimum requirements of upward and downward reserve of the system; ${\overline{Q}}_i^G$ and ${\underset{\_}{Q}}_i^G$ are, respectively, the upper and lower bounds of the reactive power output of generation units; and ${\overline{P}}_{b,t,s}^W$ and ${\overline{Q}}_{b,t,s}^W$ are, respectively, the upper limits of the active and reactive power output of the wind turbine.

(4)

Transient process constraints

The paper proposes a two-stage model that imposes constraints on both wind tur-bines and traditional generation units. In transient stability analysis, a fourth-order model is used to characterize the dynamic characteristics of the generation unit [22]. The PSASP7 governor model [23] is used as the governor and the commonly used IEEE1 excitation system model [24] is used as the excitation system. The model primarily adopts the dual-mass lumped dynamic model to describe the dynamic response of the wind turbine. This model divides the rotor and the load of the wind turbine into two mass blocks. By connecting the dynamic equations of the two mass blocks, the dynamic response of the wind turbine with wind speed and load changes can be de-scribed. The differential equations for the dual-mass lumped dynamic model are writ-ten below:

$$\frac{d\theta }{dt}={\omega }_b\left({\omega }_t-{\omega }_r\right)$$

(31)

$$\frac{d{\omega }_t}{dt}=\frac{1}{2H}\left[\frac{P_m}{{\omega }_t}-K\theta -D\left({\omega }_t-{\omega }_r\right)\right].$$

(32)

In (31) and (32), ${\omega }_t$, ${\omega }_r$, and ${\omega }_b$ are, respectively, the rotational speeds of the wind turbine, the generator, and the system synchronous speed; $\theta$ denotes the torsion angle (in radians); $K$ denotes the torsional stiffness coefficient (in p.u./rad); $D$ denotes the torsional stiffness and damping coefficient; $H$ denotes the inertia constant of the wind turbine; and $P_m$ denotes the mechanical energy converted from wind energy.

The double-fed induction generator (DFIG) model [25] can be extended into the DFIG-grid model by accounting for the grid side dynamics. If we ignore the transient behavior of rotor current and DC capacitor, the comprehensive formulation for the DFIG-grid model can be expressed as follows.

For stator circuit:

$$v_s=R_s{i}_s+\frac{d{\lambda }_s}{dt}+j{\omega }_s{\lambda }_s$$

(33)

$$p_s=v_s{i}_s^{*}+v_r{i}_r^{*}+P_{i,t,s}^W$$

(34)

$$q_s=v_{s}j{X}_{ls}{i}_s^{*}+v_{r}j{X}_{lr}{i}_r^{*}+Q_{i,t,s}^W .$$

(35)

For rotor circuit:

$$v_r=R_r{i}_r+j{\omega }_r{\lambda }_r$$

(36)

$$p_r=v_r{i}_r^{*}+v_s{i}_s^{*}$$

(37)

$$q_r=v_{r}j{X}_{lr}{i}_r^{*}+v_{s}j{X}_{ls}{i}_s^{*} .$$

(38)

For grid-side:

$$p_{i,t,s}^W=U_{j,t,s}{v}_s\mathit{sin}\left({\theta }_W-{\theta }_{grid}\right)-v_s^2{g}_{s0}-v_s{v}_r\left(g_{s1}\mathit{cos}\delta -b_{s1}\mathit{sin}\delta \right)$$

(39)

$$q_{i,t,s}^W=U_{j,t,s}{v}_s\mathit{sin}\left({\theta }_W-{\theta }_{grid}\right)-v_s^2{b}_{s0}-v_s{v}_r\left(g_{s1}\mathit{sin}\delta -b_{s1}\mathit{cos}\delta \right),$$

(40)

where $v_s$ represents the stator voltage; $v_r$ represents the rotor voltage; $i_s$ represents the stator current; $i_r$ represents the rotor current; $R_s$ and $R_r$ represent the stator resistance and rotor resistance, respectively; $X_{ls}$ and $X_{lr}$ represent the stator leakage inductance and rotor leakage inductance, respectively; ${\lambda }_s$ and ${\lambda }_r$ represent the stator flux linkage and rotor flux linkage, respectively; $p_s$ and $p_r$ represent the stator power and rotor power, respectively; $q_s$ and $q_r$ represent the stator reactive power and rotor reactive power, respectively; ${\omega }_s$ and ${\omega }_r$ represent the stator angular velocity and rotor angular velocity of wind turbine, respectively; $v_{grid}$ represents the voltage magnitude of the grid, while ${\theta }_W$ represents the voltage phase angle of the wind turbine; ${\theta }_{grid}$ represents the phase angle of the grid voltage; $g_{s0}$ and $b_{s0}$ represent the stator conductance and susceptance, respectively; $g_{s1}$ and $b_{s1}$ are the stator AC conductance and susceptance of the wind turbine, respectively; $\delta$ represents the rotor angle of the wind turbine; $P_{i,t,s}^W$ represents the active power of wind turbine; and $Q_{i,t,s}^W$ represents the reactive power of wind turbine.

(5)

Transient stability constraint

The paper analyzes the stability margin of the power system and divides the system states into stability and instability with Integrated Extended Equal Area Criterion (IEEAC). The system condition can be evaluated via the following conditions.

• Instability condition: Synchronous machine equivalent electromagnetic power and mechanical power in transient equilibrium are shown in (41). The bias of their difference is greater than 0 and tends to become larger, as shown in (42). So, it is difficult for the system to achieve dynamic equilibrium.

$$P_a\left(t_u\right)=P_{mE}\left(t_u\right)-P_{eE}\left(t_u\right)=0$$

(41)

$${\dot{P}}_a\left(t_u\right)={\left.\frac{d{P}_a}{dt}\right|}_{t=t_u}>0$$

(42)

In (41) and (42), $P_{mE}$ and $P_{eE}$, respectively, are equivalent to mechanical power and electromagnetic power; $t_u$ is the moment of reaching the unstable equilibrium; and $P_a$ is the unstable power.

• Stability condition: When the system reaches the maximum swing angle, the system’s equivalent angular velocity decelerates to 0, as shown in (44). The mechanical power of the system becomes less than the electromagnetic power, so the system will reach dynamic equilibrium, as shown in (45). The stability margin can be calculated as in (43)–(45):

$$\omega \left(t_{\overline{\theta }}\right)=0$$

(43)

$$P_t\left(t_{\overline{\theta }}\right)=P_{mE}\left(t_{\overline{\theta }}\right)-P_{eE}\left(t_{\overline{\theta }}\right)<0$$

(44)

$$f_{\sigma }=-{{\displaystyle \int }}_{\overline{\theta }}^{{\theta }_u}{P}_t\left(t_t\right)d\theta \approx \left|P_t\left(t_t\right)\right|\frac{\left({\theta }_u-\overline{\theta }\right)}{2},$$

(45)

where $t_{\overline{\theta }}$ is the moment of maximum swing angle; and $P_t\left(t_{\overline{\theta }}\right)$ and $\omega \left(t_{\overline{\theta }}\right)$ are, respectively, the unbalanced power and equivalent velocity of the maximum swing angle.

The system stability information under the fault scenario can be obtained by using Monte Carlo sampling. If the system meets the constraint (43)–(45), it is in a stable condition, namely $n_{ss}=1$; if the system satisfies the constraint (41) and (42), the system is in an instable condition. Through the analysis of S samples, the probability information of system transient stability can be obtained.

3.2. Emergency Control Stage Model

3.2.1. Objective Function

In case of a system fault, it is advisable to apply Formulae (41)–(45) to ascertain whether the system has lost its transient stability. If the system is operating steadily in line with (43)–(45), then there is no urgent need for emergency control measures. However, if the system is in an unstable status that is in line with (41) and (42), emergency control measures become necessary. In this case, fast response of emergency control measures is extremely important.

In the emergency control stage, the optimization is carried out with the goal of minimizing the operational costs of the system. The objective is shown in (46), which contains two parts: the operation cost of generation units (47) and the operation cost of energy storage (48), both of which include the cost of delivering reserve. The uncertainty in the emergency control stage makes the problem difficult to solve; thus, in the emergency control stage, only the forecasted scenarios are considered. After the fault is cleared, the system returns to the preventive control stage using stochastic analysis.

$$\min f_3=f_{G1}+f_{S1}$$

(46)

$$f_{G1}={{\displaystyle \sum }}_{i=1}^{n_G}[r_{i,t}^{G,UR}{C}^{UR}+r_{i,t}^{G,DR}{C}^{DR}]$$

(47)

$$f_{S1}={{\displaystyle \sum }}_{i=1}^{n_S}\left(r_{i,t}^{RC}{C}_{i,t}^{RC}+r_{i,t}^{RD}{C}_{i,t}^{RD}\right)$$

(48)

In (46) and (47), $r_{i,t}^{G,UR}$ and $r_{i,t}^{G,DR}$ are, respectively, the delivered upward and downward reserve from generation units; $C^{UR}$ and $C^{DR}$ are the costs to deliver the upward reserve and downward reserve, respectively; similarly, $r_{i,t}^{RD}$ and $r_{i,t}^{RC}$are, respectively, the delivered upward and downward reserve from energy storage systems; and $C_{i,t}^{RD}$and $C_{i,t}^{RC}$ are the costs to deliver the upward reserve and downward reserve, respectively.

3.2.2. Constraints

During the emergency control stage, the system still needs to satisfy the operational constraints mentioned in the preventive control stage. Only by satisfying the above constraints can the system ensure stable operation and optimal economy after the disturbance.

Unlike the preventive control stage, the power balance in the emergency control stage considers the delivered reserve from generation units and energy storage; thus, it is rewritten as in (49) and (50).

$$p_{b,t}^G+p_{b,t}^W-p_{b,t}^C+p_{b,t}^D-p_{b,t}^L+r_{b,t}^{G,UR}-r_{b,t}^{G,DR}+r_{b,t}^{RD}-r_{b,t}^{RC}-u_{b,t}{{\displaystyle \sum }}_{j=1}^{n_b}{u}_{j,t}\left(G_{bj}cos{\delta }_{bj,t}+B_{bj}sin{\delta }_{bj,t}\right)=0$$

(49)

$$q_{b,t}^G+q_{b,t}^W-q_{b,t}^L+u_{b,t}{{\displaystyle \sum }}_{j=1}^n{u}_{j,t}\left(B_{bj}cos{\delta }_{bj,t}-G_{bj}sin{\delta }_{bj,t}\right)=0$$

(50)

Different from the preventive control stage, reserved generation capacity during the emergency control stage will be effectively delivered. Thus, the power limit constraint should be revised as in (51). Since the upward and downward reserve capacities have been determined in the preventive control stage, the available reserve capacities in the emergency control stage cannot exceed the determined amounts. These limits are enforced by (52) and (53). The constraints concerning the reactive power outputs of generator units and the active/reactive power outputs of wind turbines remain the same as shown in (54)–(56). The node voltage and the transmission line power flow constraints, as shown in (57) and (58), should be strictly satisfied.

$${\underset{\_}{P}}_i^G+r_{i,t}^{G,DR}\le p_{i,t}^G\le {\overline{P}}_i^G-r_{i,t}^{G,UR}$$

(51)

$$r_{i,t}^{G,UR}\le {\mathcal{R}}_{i,t}^{G,UP}, r_{e,t}^{RC}\le {\mathcal{R}}_{e,t}^C$$

(52)

$$r_{i,t}^{G,DR}\le {\mathcal{R}}_{i,t}^{G,DN}, r_{e,t}^{RD}\le {\mathcal{R}}_{e,t}^D$$

(53)

$${\underset{\_}{Q}}_i^G\le q_{i,t}^G\le {\overline{Q}}_i^G$$

(54)

$$0\le p_{b,t}^W\le {\overline{P}}_{b,t}^W$$

(55)

$$0\le q_{b,t}^W\le {\overline{Q}}_{b,t}^W$$

(56)

$${\underset{\_}{U}}_b\le u_{b,t}\le {\overline{U}}_b$$

(57)

$$-{\overline{S}}_l\le s_{l,t}\le {\overline{S}}_l$$

(58)

The constraints of energy storage systems in the emergency control stage remain the same as in the preventive control stage, as shown in (59) and (60). The reserve capacity constraint for energy storage is shown in (61). Similar to the generation units, the power and stored energy of the energy storage will change after delivering the re-served energy in the emergency control stage. Its modified model is formulated as in (59)–(63). It is worth noting that the stored energy evolution constraint has been adjusted considering the delivery of reserved energy. The energy evolution constraint is revised as in (61).

$${\underset{\_}{P}}_e^C{c}_{e,t}\le p_{e,t}^C+r_{e,t}^C\le {\overline{P}}_e^C{c}_{e,t}$$

(59)

$${\underset{\_}{P}}_e^D{d}_{e,t}\le p_{e,t}^D+r_{e,t}^D\le {\overline{P}}_e^D{d}_{e,t}$$

(60)

$$e_{e,t}^S=e_{e,t-1}^S+\left(p_{e,t}^C+r_{e,t}^C\right){\eta }_e^C\Delta t-(p_{e,t}^D+r_{e,t}^D)/{\eta }_e^D\Delta t$$

(61)

$${\underset{\_}{E}}_e^S\le e_{e,t}^S\le {\overline{E}}_e^S$$

(62)

$$c_{e,t}+d_{e,t}\le 1$$

(63)

4. Solution Method

The two-stage model proposed in this paper involves two different types of problem, and therefore requires different algorithms. The preventive control model is a multi-objective collaborative optimization problem, which is handled by a multi-objective particle swarm optimization (MPPSO) algorithm; and the model of the emergency control model is a single-objective problem, which is solved by a particle swarm optimization algorithm (PSO).

The multi-objective particle swarm algorithm solves the multi-objective problem with Pareto front solutions [26], which is illustrated in Figure 1.

4.1. Initialization

The single-objective PSO algorithm and the MOPSO algorithm have the same rules in particle initialization and iteration. Usually, the initialization of particle swarm is conducted to randomly generate a group of particles as the initial solution. After randomly generating the initial particles, if they are in the inferior solution space, the optimization results may be compromised. In this paper, when solving the advanced preventive control model, the operation schedule data of the previous day are used to set the initial values of particles; and when solving the emergency control model, the current operation status is used for the initial values of particles.

4.2. Particle Update

The location and velocity of PSO and MOPSO algorithms can be calculated as in (64):

$$\left\{\begin{array}{c}{x}_i\left(t+1\right)=x_i\left(t\right)+v_i\left(t+1\right)\\ v_i\left(t+1\right)=\omega \left(t\right)v_i\left(t\right)+c_1{r}_1\left(pbest-x_i\left(t\right)\right)+c_2{r}_2\left(gbest-x_i\left(t\right)\right)\end{array}\right.$$

(64)

where $x_i\left(t+1\right)$ and $v_i\left(t+1\right)$ are, respectively, the location and velocity of particle $i$ in the ${\left(t+1\right)}^{th}$ iteration; $gbest$ is the global best position of the particle; and $\omega$ is the inertial factor whose magnitude determines the strength of the global and local search ability. In (64), $v_i$ is the memory term, which denotes the velocity of particle $i$; $x_i$ is the current position of particle $i$; $c_1$ and $c_2$ are learning factors and, usually, $c_1=c_2=2$; $rand$ denotes the random number between (0, 1); $r_1$ represents the individual cognitive factor, which is used to adjust the influence of a particle’s personal best position during the update; $r_2$ represents the social cognitive factor, which is used to adjust the influence of the global best position during the update; $pbest$ represents the personal best position of a particle; and $gbest$ represents the global best position—it represents the best solution found by all particles.

To accelerate the convergence of the algorithm, a method that decreases the inertial velocity of particles, called linearly decreasing weight (LDW), is used as in (65).

$$\omega \left(t\right)=\frac{\left({\omega }_{ini}-{\omega }_{end}\right)\left(G_k-t\right)}{G_k}+{\omega }_{end}$$

(65)

In (65), $G_k$ denotes the maximum number of iterations; ${\omega }_{ini}$ denotes the initial inertia weight, which is typically set as 0.9; and ${\omega }_{end}$ denotes the inertia weight when the iteration number reaches the maximum, which is typically set as 0.4.

4.3. Fitness Function and Optimal Solution Construction

For the advanced prevention control model, the weights of the two different objectives need to be determined. By contrast, the emergency control model is a single-objective optimization problem, so Pareto optimal solution is not involved. In this paper, the entropy method [27] and TOPSIS [28] method are used to deal with the weight of multiple objectives. The combination of the two can help effectively find the objective weight of the Pareto optimal frontier, and then the TOPSIS method can sort these Pareto optimal solutions as options. The steps are as follows:

(1)

Obtain the original indicator data matrix based on sample data. Let a matrix with n-dimensional objective function and$m$ Pareto solutions be $X={\left(x_{ij}\right)}_{m\times n}$; then, standardize it to obtain a new matrix $Y={\left(y_{ij}\right)}_{m\times n}$:

$$\left\{\begin{array}{c}{y}_{ij}=\frac{x_{ij}-\mathit{min}{x}_j}{\mathit{max}{x}_j-\mathit{min}{x}_j} gain index\\ y_{ij}=\frac{\mathit{min}{x}_j-x_{ij}}{\mathit{max}{x}_j-\mathit{min}{x}_j} cost index\end{array}\right.$$

(66)

where $\mathit{max}{x}_j$ is the maximum value of all evaluation data of index $j$; and $\mathit{min}{x}_j$ is the minimum value of all evaluation data of index $j$. The gain index indicates that the larger the value, the better the index. The cost index is the opposite.

(2)

Obtain the evaluation indicator $p_{ij}$ and entropy $E_j$ of samples; then, obtain the information redundancy value $d_j$ based on the information entropy, as follows:

$$p_{ij}=\frac{y_{ij}}{{{\displaystyle \sum }}_{i=1}^m{y}_{ij}}$$

(67)

$$E_j=-\ln n^{-1}{\displaystyle \sum }_{i=1}^m{p}_{ij}\ln p_{ij}$$

(68)

$$d_j=1-E_j .$$

(69)

(3)

Finally, calculate the weight of each index according to the information entropy:

$$w_j=\frac{d_j}{{{\displaystyle \sum }}_{j=1}^n{d}_j} .$$

(70)

(4)

Standardize the decision matrix to obtain the normalized attribute value of the original scheme and introduce the weight vector $V^{+}=\left(V_1^{+},V_2^{+},\dots V_n^{+}\right)$ obtained via the entropy weight method:

$$z_{ij}=\frac{x_{ij}}{\sqrt{{{\displaystyle \sum }}_{i=1}^m{x}_{ij}{}^2}}$$

(71)

$$V_{ij}=w_j{z}_{ij} .$$

(72)

(5)

Determine the optimal solution $V^{+}$ and worst solution $V^{-}$:

$$\left\{\begin{array}{c}{V}^{+}=\left(V_1^{+},V_2^{+},\dots V_n^{+}\right)\\ V^{-}=\left(V_1^{-},V_2^{-},\dots V_n^{-}\right)\end{array}\right.$$

(73)

(6)

Calculate the distance from each scheme to the optimal solution and the worst solution in (74) and (75); then, calculate the relative proximity of each scheme to the ideal solution according to Formula (75):

$$\left\{\begin{array}{c} D_i^{-}=\sqrt{{\displaystyle \sum }_{j=1}^n{\left(V_{ij}-V_j^{-}\right)}^2}\\ D_i^{+}=\sqrt{{\displaystyle \sum }_{j=1}^n{\left(V_{ij}-V_j^{+}\right)}^2}\end{array}\right.$$

(74)

$$O_i=\frac{D_i^{-}}{D_i^{+}+D_i^{-}}$$

(75)

The flowchart of the entropy-TOPSIS method is shown in Figure 2.

5. Case Study

5.1. System Discerption

The proposed model is validated with the modified IEEE 118-node testing system, which is shown in the Figure A1, in Appendix A. The system data are available at [29]. The modified IEEE 118-node system contains ten wind farms, each of which contains 100 double-fed wind turbines (2 MW). The rated power and capacity of the energy storage systems are, respectively, 10 MW and 20 MWh. The costs of upward and downward reserve capacities are set as CNY 20/MW. The locations of the wind farms and energy storage systems are shown in Appendix A. The system voltage limits are set as [0.97, 1.06] p.u. The uncertainty of the system is simulated via Monte Carlo simulation. The corresponding distributions and parameters are set as follows:

(1)

The nodal load power is assumed to follow a normal distribution with a mean value given by the provided standard data and a variance of ±10%. The load correlation coefficient of the connected node is set to 0.15.

(2)

The cut-in, cut-out, and rated wind speeds of wind turbines are, respectively, set as 3 m/s, 25 m/s, and 12 m/s.

(3)

To simulate a fault on branches 18–35, the probability values of fault type and location are adopted from Table 1 and

Table 2. Costs of preventive control stage.

Cases	Case 1	Case 2	Case 3
Costs (USD/h)	95,634	98,985	104,182

. The fault recovery time is assumed to follow a normal distribution with a mean value of 350 ms and a standard deviation of 10%.

The non-sequential Monte Carlo sampling method is used to generate 1000 scenarios, and K-MEANS clustering is used to reduce them to 15 typical scenarios.

This paper presents a rolling simulation analysis using the provided case study distribution starting at 10:00 a.m. During the simulation, generation unit 17 (40 MW) experiences a fault at 12:00 p.m. and recovers to the stability condition afterward. At 1:00 p.m., generation units 2, 10, 16, 21, 27, 35, and 43 experience faults. Specifically, the system is in the preventive control stage between 10:00 a.m. and 12:00 p.m. Although a minor fault occurs at 12:00 p.m., the system remains stable after fault recovery, maintaining preventive control. After 2:00 p.m., the system enters the emergency control stage. The time steps in the preventive and emergency control stages are, respectively, 15 min and 0.1 s. The solving times of the two states are 2 h and 8 s. In addition, to further analyze the impacts of the energy storage systems and the reserve requirements of the system, three different system configuration schemes are used for analysis:

• Case 1: Reserve and energy storage systems are not considered;
• Case 2: Only reserve from generation units is considered (energy storage systems are not considered);
• Case 3: Reserve from both generation units and energy storage systems are considered.

5.2. Optimization Results

5.2.1. Analysis before Failure

From 10:00 a.m. to 11:59 a.m., the system remains fault-free and is in the preventive control stage. In Case 1, the maintenance work of the system remained in the non-optimized state as in the previous moment, and it was still undergoing the preventive control stage. The average cost per hour in this scenario is USD 95,634/h.

In Case 2, although the preventive control stage is initiated, the energy storage participation is not considered, and only the reserve capacity from generation units is considered. Due to the inclusion of the reserve cost of the generation units, the total cost before the fault increases to USD 98,985/h.

In Case 3, the reserve from energy storage systems is considered in the preventive control stage, which causes an increase in the cost before the fault to USD 104,182/h.

In the preventive control stage, the scheduled reserve capacities for Case 2 and Case 3 are illustrated in Figure 3 and Figure 4, respectively.

During the preventive control stage, the system’s demand for the reserve is mainly concentrated on the peak load hours between 11:00 and 12:00. Based on Figure 4 and Figure 5, it can be observed that the scheduled reserve capacities in Case 2 and Case 3 in the preventive control phase can meet the reserve capacity demand during peak load hours. In Case 3, the energy storage systems provide 4.89 MW of upward reserve capacity and 3 MW of downward reserve capacity in the loop. A horizontal comparison reveals that the energy storage systems can provide a certain amount of upward re-serve capacity during peak load hours, fulfilling about 5.4% of the system’s reserve capacity requirement. This would enable other power supply units to remain at high output levels without increasing the operational risk of the system. During valley load hours, the energy storage system can provide a certain amount of downward reserve capacity to ensure the optimal utilization of resources.

The probabilities of the system stability with the three cases are shown in Figure 5.

In this scenario, the other probability indicators of the system do not meet the requirements; for instance, in the case of a three-phase short circuit fault with a fault lo-cation near node 22, ① the probability of node 21 voltage magnitude being in the range of [0.97, 1.06] is only 0.2, which is lower than the required probability of 0.95; ② the probability of the transmission line flows being within the limit is only 0.6; and ③ the probability of the system attaining transient stability is 0.54. Additionally, there are differences in safety margins among the three schemes. In Case 1, for probability indicators, such as the probability of node voltage magnitude, the probability of the trans-mission line flows being within the limit and the probability of then attaining transient stability are both below the threshold. As a result, the obtained safety margin is only 0.52. In Case 2 and Case 3, the probability indicators of the system are improved to a certain extent, and their stability margins reach 0.95 and 1.00, respectively. Particularly in Case 3, by optimally controlling the reserve capacity and considering energy storage systems’ participation, the system is better secured.

5.2.2. Analysis on Single Failure

At 1:00 p.m., generation unit 17 (40 MW) encounters a fault. The power angle trajectories of the generation unit in all three cases are shown in Figure 6.

In Case 1 and Case 2, the system loses stability in 2–4 s after the fault occurs. Case 1 loses stability first due to the lack of reserve support. In contrast, the re-serve capacity of Case 3 was adequate to make up for the gap caused by the fault. The power angle oscillated within 50 degrees, and the power angle deviation gradually de-creases over time. The system shows a tendency towards stability recovery.

5.2.3. Analysis on Multiple Faults

At 2:00 p.m., generation units 2, 10, 16, 21, 27, 35, and 43 experience faults. With-out implementing emergency control measures, the rotor angle trajectories of the generation units in all three cases are shown in Figure 7.

As shown in Figure 7, in all three cases, the system becomes unstable after the faults occur in the absence of emergency control measures. Like the single fault case, the system in Case 1 becomes unstable first; by contrast, the system in Case 3 can re-main stable for some time, relying on the reserve from the generation units and the energy storage systems, but it will eventually become unstable. The rotor angle trajectories for generation units subjected to emergency control measures are shown in Figure 8.

The system in Case 1 will still become unstable after 4 s of emergency control as it does not account for the reserve capacity from generation units and energy storage systems, which results in significant power shortages during multiple faults. In such a case, the system cannot compensate for the power shortage; hence, it loses stability. By contrast, Case 2 considers the reserve capacity from generation units, which helps reduce the deviation of the generator power angle within 8 s after the emergency control, indicating a trend towards stable operation. In comparison, for Case 3, which further considers the reserve from energy storage system, the deviation of the generator power angle is reduced by approximately 10%, and it exhibits a quicker recovery to stable operation. In Appendix A, Figure A2 and Figure A3 illustrate the reserve adjustment schemes.

The costs of the three cases in the emergency control stage are shown in

Table 3. Costs of the emergency control stage of the three cases.

Cases	Case 1	Case 2	Case 3
Costs ($/h)	98,421	92,112	90,101

.

Among the three cases, Case 1 incurs the highest cost since it fails to account for the reserve capacity of generation units and the energy storage system, resulting in difficulty stabilizing the system. In contrast, Case 2 considers the cost of reserve during the preventive control stage, and the inclusion of reserve from generation units helps the system to remain stable, which reduces the cost by 5.1% compared to that of Case 1. Moreover, Case 3 further considers the reserve capacity from energy storage systems and has the lowest adjustment cost; it reduces the cost by 9.6% compared to Case 1 and by 2.2% compared to Case 2.

5.2.4. Algorithm Comparison

To demonstrate the effectiveness of the proposed improvement to the solving algorithm, Case 3 is taken as an example. The improved MOPSO (IMOPSO) and the traditional MOPSO algorithms are used to solve the optimization problem for each objective. The results are compared in

Table 4. Results of different optimization algorithms.

Algorithm	Costs ($)	Iterations	Average Calculation Time for Each Iteration
MOPSO	98,613	183	4.52 s
IMOPSO	90,101	145	4.43 s

.

Table 4 reveals that the comprehensive cost obtained by the IMOPSO algorithm is 8.63% lower than that obtained by the traditional MOPSO, while the iteration speed of IMOPSO is 20.76% fast than that of the traditional MOPSO, demonstrating the effectiveness of the algorithm improvement proposed in this paper, which enhances the global search capability of the algorithm.

6. Conclusions

This paper proposes a two-stage preventive control optimization model that ad-dresses the possible system transition instability issue caused by the uncertainties introduced by different system components considering the balance of economy and system safe operation. The first stage of the model is used to obtain an optimal operation schedule which guarantees a high probability of the system being stable and a low operation cost. The second stage is applied to respond to the emergency after the occurrence of any faults that may result a stability issue. Furthermore, the multi-objective particle swarm optimization algorithm based on entropy-TOPSIS is adopted to accelerate the solving process of the two-stage model. The case study shows that with the proposed model, while pursuing economics, the safety and stability mar-gin of the system can be reduced. In addition, when multiple faults occur and threaten the system’s stability, emergency control can help the system meet the satisfy the stability margin requirements with a slightly higher cost on the top of the schedule obtained in the preventive control stage.