Integrated Transmission Network Planning by Considering Wind Power’s Uncertainty and Disasters

Abstract

The penetration of wind turbines and other power sources with strong uncertainty into the grid has increased in recent years. It has brought significant technical challenges to power systems’ operation. The volatility and intermittency of wind power increase the risk of insufficient transmission capacity of the lines. Therefore, the traditional deterministic planning methods for transmission grids are no longer fully applicable. On the other hand, the frequent disasters in recent years have posed a great threat to the power system, especially for the transmission grid. This requires the design of transmission lines with high design standards, such as skeleton networks, to withstand disasters. With the aim to address these problems, a bi-level integrated network planning model for the transmission grid is developed by considering wind power’s uncertainty and load guarantee under disasters. Chance constraints are used in the model to characterize wind power’s uncertainty, and a skeleton network is adopted to cope with disasters. Moreover, based on a convex relaxation method, the chance constraints are converted into the probabilistic inequalities to be solved. The proposed method is simulated in the IEEE 118 bus system, and the obtained network planning scheme is further analyzed in the scenario tests. And the result of the tests proves the validity and reasonableness of the proposed method.

1. Introduction

In response to the growing demand for electricity, the power system has entered a new stage of developing large grids and large units [1]. On the other hand, renewable energy generation, represented by wind power, can effectively reduce the emission of greenhouse gases. This brings a non-negligible contribution to global warming. Consequently, the increased penetration of renewable energy sources is the general trend. However, the high degree of uncertainty is a distinctive feature that distinguishes it from conventional generators [2]. All these have put forward higher requirements in the network planning. As a key component of the transmission planning, network planning is significant for the safe and reliable operation of the power system [3].

It is worth noting that the frequency of extreme disasters has increased more in recent decades. This has caused greater hazards to infrastructures with more exposed equipment, such as transmission grid [4]. For example, the extreme ice disaster that occurred in South China in 2008 caused more than 450 lines to be blacked out. It resulted in a total of more than 14.66 million residents suffering from power outages [5]. This indicates that the previous standard of lines has not been able to cope with the increasingly frequent extreme disasters. Therefore, the network for transmission needs to be strengthened to withstand extreme disasters in the network planning [6].

The network planning for transmission refers to the design of grid structure to meet the requirements of safe operation within a given planning period [7]. The basic principle of network planning is to ensure the reliability and stability of the power system, while making the network construction and operation costs minimal [8]. With the increasing penetration of renewable energy, some renewable energy sources cannot be used effectively due to the weakness of the grid structure. This leads to the forced abandonment of this part of energy. Taking wind power as an example, wind is highly random and intermittent, so the wind power is also highly uncertain. Furthermore, the lines are not able to withstand the risk of power flow crossing the limits, which are caused by the uncertainty of wind power. In some cases, even the scheduling optimization of the turbines cannot improve these situations [9].

Considering the more extreme case, some lines would stop operation if serious disasters occur. However, the transmission network must guarantee the supply of critical loads in the system through the dangerous period of regional blackouts [10]. Improving the disaster prevention standard of lines is a way to enhance the defense performance for the grid. However, different lines of the transmission are required to carry different load profiles and are in different positions in the transmission process. Planning with enhancement of all lines would be more of an investment and less economical [11]. Thus, the skeleton network for the transmission grid needs to be established. The skeleton network is relied upon to secure the supply of critical loads and to lay the foundation for subsequent system maintenance to restore power [12].

Many studies have been carried out on the network planning for the transmission. The network planning problems are essentially mixed integer, nonlinear, nonconvex optimization problems [13]. In general, when the system’s size increases, the model becomes more complex and requires a large computational memory. Therefore, researchers usually take some methods to reduce the size of the problem or to get the fast optimal solutions:

• Heuristic algorithms [14] or metaheuristic strategies [15,16] are used for solving the problem. But it has problems such as high processing requirements and difficulty in specifying the stopping criteria.

• Use modeling methods such as hybrid models [17] and parsimonious models [18] to reduce the complexity brought by the nonlinearities, but there is a risk of running into trouble with local optimal solutions.

• Use intelligent algorithms such as genetic algorithms [19]. This type of method can correspond well to the practical significance of line construction. However, for larger scale systems, the computation time grows exponentially.

Existing studies have incorporated the examination of wind power’s uncertainty while conducting transmission planning. The current main methods can be divided into two categories: scenario analysis [1,20,21,22,23,24] and mathematical analysis [25]. The literature [20] uses scenario reduction techniques for wind power’s uncertainty. Although the scenario reduction can effectively reduce the size of the problem, the selected scenario greatly determines their accuracy. In the literature [25], the probability density function of wind power is considered as the DC model based on the chance constrained planning. And the optimal solution is found under the confidence level, but the calculation of its convolution is complicated, and the fitting effect of wind power needs to be improved.

In general, most of the existing studies have unilaterally considered the impact of wind power’s uncertainty or disasters on transmission network planning, but in fact the two are related. Therefore, in the planning process of the grid, the impact of wind power should be considered on the one hand, and the original or new lines should be enhanced differently to improve the disaster resistance on the other hand. Therefore, an integrated transmission network planning method is proposed, which takes wind power’s uncertainty and grid resilience into account together, with the following contributions:

• A bi-level network planning model is established to divide the problem into a decision level and an operation level, which effectively reduces the model’s complexity.

• The uncertainty of wind power is expressed in the form of chance constraint. They are converted into the probability inequalities by the convex relaxation method. That is useful for the solution under different confidence probabilities.

• The “big-M” method is used to separate the 0–1 decision variables, which effectively deals with the mixed-integer planning problem of integrated network planning.

The rest of the paper is organized as follows: Section 2 introduces the fitting method of wind power’s uncertainty with Gaussian mixture model. Section 3 introduces the responsibilities of the skeleton network and the framework for integrated network planning. Section 4 develops a bi-level model of integrated network planning for the transmission. Section 5 transforms the chance constrained planning model based on the convex relaxation method. In Section 6, the effectiveness of the proposed method is verified by simulation tests on the IEEE 118 bus system. The conclusion is drawn in Section 7.

2. Wind Power’s Uncertainty and Gaussian Mixture Model

Wind turbines are affected by climate, geography, and other factors, and their power contains fluctuation and intermittency. Therefore, the probability distribution of wind power has strong uncertainty, which is difficult to be described uniformly by some analytical function. Most studies use single probability distributions such as the Normal distribution, Beta distribution [26], and so on to approximate the uncertainty. However, it will cause statistical errors. In contrast, the Gaussian mixture model (GMM) can more accurately fit the irregular or multi-peaked characteristics of wind power. And it can accommodate more complex and diverse distributions. Therefore, GMM has more applications in dealing with uncertainties and other problems [26,27].

2.1. Gaussian Mixture Model

GMM is based on a weighted convex combination of multiple Gaussian distribution functions [27] to describe the probability density distribution of random variables. For a wind power plant j, its wind power output ${\tilde{P}}_j^W$ is considered as the sum of the output forecast ${\overline{P}}_j^W$ and the output forecast deviation ${\tilde{\epsilon }}_j^W$, that is ${\tilde{P}}_j^W={\overline{P}}_j^W+{\tilde{\epsilon }}_j^W$. ${\overline{P}}_j^W$ is obtained by the regulation center based on a large amount of historical data of wind power and meteorological conditions observed before dispatch. While ${\tilde{\epsilon }}_j^W$ is a random variable. GMM can be classified as

(1)$H({\tilde{\epsilon }}_j^W)={\displaystyle \sum _{m=1}^M{\omega }_{m,j}}{\varphi }_{m,j}\left({\tilde{\epsilon }}_j^W\mid {\mu }_{m,j},{\sigma }_{m,j}^2\right)$

where M is the number of Gaussian components; ω_m,j is the weight of the m_th Gaussian component and ${\displaystyle \sum _m{\omega }_{m,j}=1}$; ϕ_m,j is the m_th Gaussian component of the output of wind power plant j; μ_m,j and ${\sigma }_{m,j}^2$ are the expectation and variance corresponding to this Gaussian component, respectively.

For the aim to construct a GMM with wind power, it is crucial to obtain ω_m,j, μ_m,j, and ${\sigma }_{m,j}^2$ of each Gaussian component. The generally accepted Expectation-Maximum (EM) algorithm is an iterative strategy for parameter estimation for the above parameters. Each iteration is divided into two steps: E-step and M-step.

E-step: Calculate the probability or weight ω_im of a data point from the m_th Gaussian component from a large amount of data based on the current parameters.

(2)${\omega }_{im}=\frac{{\omega }_m{\varphi }_m({\epsilon }_i^W|{\mu }_m,{\sigma }_m^2)}{{\displaystyle {\sum }_{r=1}^M{\omega }_r{\varphi }_m({\epsilon }_i^W|{\mu }_m,{\sigma }_m^2)}}$

M-step: Parameter estimation for this round of iteration is based on the E-step:

(3)${\omega }_m=\frac{{\displaystyle {\sum }_{i=1}^N{\omega }_{im}}}{N}$

(4)${\mu }_m=\frac{{\displaystyle {\sum }_{i=1}^N{\omega }_{im}{\epsilon }_i^W}}{{\displaystyle {\sum }_{i=1}^N{\omega }_{im}}}$

(5)${\sigma }_m^2=\frac{{\displaystyle {\sum }_{i=1}^N{\omega }_{im}{({\epsilon }_i^W-{\mu }_m)}^2}}{N}$

In Equations (2)–(5), N is the total number of data; ${\epsilon }_i^W$ is the i_th data. The obtained parameter estimates can be considered as valid feasible solutions when the number of iterations for the EM algorithm is large enough [27].

2.2. AIC Evaluation

The advantage of GMM is that it is not limited to a specific probability distribution. Theoretically, it can fit any type of probability distribution with high accuracy when M of GMM is large enough. However, if M is too large, the complexity of the model increases and even leads to overfitting. Therefore, the Akaike Information Criterion (AIC) [28] is introduced to effectively evaluate and determine the suitable M:

(6)$\mathrm{AIC}=2k-2\ln (L({\epsilon }^W))$

where k is the total number of parameters to be estimated; $L({\epsilon }^W)$ is the similarity function of GMM. Therefore, the AIC can be rewritten as:

(7)$\mathrm{AIC}=2k-2{\displaystyle {\sum }_{i=1}^N\ln ({\displaystyle {\sum }_{m=1}^M{\omega }_m{\varphi }_m({\epsilon }_i^W|{\mu }_m,{\sigma }_m^2))}}$

The reasonable GMM of the wind power is constructed by selecting the appropriate M by minimizing AIC according to the above equations.

3. Skeleton Network’s Responsibilities and Framework of Integrated Network Planning

3.1. The Responsibilities of Skeleton Network against Disasters

The integrated network planning for the transmission should not only enable the system to effectively accommodate a series of changes caused by the uncertain renewable energy under normal operating. It should be also resilient to the extreme disasters.

The inherent structure of the transmission and its geographical location make it vulnerable to disasters. If the lines connecting to the critical loads are out of operation due to disasters, the supply of critical loads may be difficult to secure. Selecting and constructing a skeleton network from the transmission is an effective way to increase the resilience against disasters [10]. A skeleton network against disasters is a smaller core backbone (e.g., the red bolded lines in Figure 1) that meets topological connectivity and safe operation constraints by increasing the design criteria of some lines. The skeleton network can be used to secure critical loads under extreme disasters. The critical loads guaranteed are the loads necessary to ensure social, political, and economic stability. The critical loads generally cover emergency response load for society such as hospital, firefighting, transportation, electricity repair, and basic residential electricity; critical pivotal load such as communication equipment, banks, and government offices; and high-risk customer loads such as the toxic chemical industry, mineral industry, and nuclear industry [7].

In addition, considering the huge investment cost is required to enhance the construction of all lines without any difference. Therefore, it is necessary to screen out important lines for differential enhanced planning according to the critical loads. And these lines can form a skeleton network to guarantee the supply of critical loads.

3.2. The Framework of Integrated Network Planning for the Transmission

As the penetration of fluctuating renewable power such as wind power continues to increase, the transmission grid’s ability to withstand fluctuating power also needs to be increased. Currently, there are many measures to deal with the power flow overload of the lines by dispatching generators or abandoning wind power. However, the scope of the above regulations is still limited. And the new lines should be built to further enhance the ability to withstand the uncertainties brought by wind power.

And when extreme disasters occur, the skeleton network against disasters should satisfy the reliable and continuous power supply for the critical loads. At the same time, the normal loads should also be satisfied as much as possible.

Figure 2 is a schematic diagram of the framework of integrated network planning for the transmission. New lines need to be built to accept wind power under normal operation. And a collection of lines to be enhanced is selected to form a skeleton network to cope with the supply of critical loads under disasters (the red nodes indicate that they are connected to critical loads). The lines of the skeleton network are enhanced based on the ordinary lines. And after enhanced construction, the power flow capacity of these lines will also be further increased, thus affecting the planning results of under normal operation. Therefore, the new ordinary line planning and the skeleton network planning are mutually influential, and the two together constitute a comprehensive network planning for the transmission.

4. A Bi-Level Model for Integrated Network Planning for the Transmission

4.1. Bi-Level Planning Model

The transmission network planning problem is a mixed integer nonlinear problem [13]. It is difficult to solve directly, especially for larger systems. The model is more complex when uncertainties such as wind power are considered. Therefore, suitable planning model structures need to be studied to reduce the solving difficulty. The bi-level planning model is an effective method to solve such complex problems, and it has been used in many articles of research for power systems [17].

The general form of the bi-level planning model is:

(8)$\min F=F(x,v)$

(9)$s.t.G(x)\le 0$

(10)$\min v=f(x,y)$

(11)$s.t. \mathrm{g}(x,y)\le 0$

where F(x,v) and $G(x)\le 0$ are the objective function and the constraints of upper-level model, respectively; f(x,v) and $\mathrm{g}(x,y)\le 0$ are the objective function and the constraints of lower-level model, respectively; v is the objective value of lower-level model; x and y are the decision vector of upper and lower-level model, respectively.

In the bi-level planning model, the upper-level decisions influence the lower-level objectives and constraints, which affect the lower-level decisions. The lower-level in turn returns with the optimal objective value to influence the upper-level decisions. This realizes the interaction between the upper and lower levels.

4.2. Integrated Network Planning Model for the Transmission

4.2.1. Upper-Level Model

The bi-level planning model is used to integrate the network planning problem by layering to make comprehensive decisions. The upper-level model is the network structure scheme layer, and its objective function is Equation (12). It mainly includes five terms: the investment cost of new lines, the cost of enhancing ordinary lines, the penalty for the predicted wind power’s abandonment, the supply benefit of critical loads, and the supplementary penalty of wind power’s abandonment returned from the lower-level model.

(12)$\min F={\displaystyle \sum _{t\in S_{plan}}[C_{nor}{H}_t{L}_t}+C_{bone}{D}_t{L}_t]+{\displaystyle \sum _{t\in S_{ori}}{C}_{bone}{D}_t{L}_t}+{\displaystyle \sum _{j\in {\varphi }_w}\alpha \Delta p_j^W}-{\chi }_L{\displaystyle \sum _{i\in S_N}{P}_{L_i}^{(d)}}+\delta f(\Delta p_{ad}^W)$

where S_plan is the set of optional construction lines; S_ori is the set of original lines; ϕ_w is the set of buses connected with wind power; S_N is the set of all buses of the system; C_nor and C_bone are the unit costs of building lines as ordinary lines and enhanced lines, respectively; H_t, D_t represent the 0–1 decision variables of whether to build line t as ordinary lines or enhanced lines, respectively. A value of 0 means no construction, and 1 means construction; L_t is the length of line t; $\Delta p_j^W$ and α are the predicted wind power’s abandonment and its penalty coefficient, respectively; $P_{L_i}^{(d)}$ is the loads to maintain supply, and the superscript (d) indicates that it is under disasters; χ_L is the per benefit for supplying power. The more important the load is, the larger its value is. The guaranteed load supply benefits should be taken to maximize, so this item is taken to be negative; $f(\Delta p_{ad}^W)$ is the optimization result returned by the lower-level model; δ is its penalty coefficient.

The upper-level model is subject to the following constraints for the normal operation (denoted by superscript (0)) and the state under disaster (denoted by superscript (d)):

s.t.

(13) $P_{Gi}^{(0)}+{\overline{P}}_i^W-\Delta p_i^W-P_{L_i}^{(0)}-{\displaystyle \sum _{ij\in S_{ori}}{p}_{ij}^{(0)}}-{\displaystyle \sum _{ij\in S_{plan}}{\widehat{p}}_{ij}^{(0)}}=0,\forall i\in S_N$

(14) $p_{ij}^{\left(0\right)}-({\theta }_i^{\left(0\right)}-{\theta }_j^{\left(0\right)})/x_{ij}=0,\hspace{1em}\forall ij\in S_{ori}$

(15) $\{\begin{array}{c}{p}_{ij}^{\left(0\right)}\ge (1-{\lambda }_l)[-P_{nor}^{\max }-D_{ij}(P_{bone}^{\max }-P_{nor}^{\max })]\\ p_{ij}^{\left(0\right)}\le (1-{\lambda }_l)[P_{nor}^{\max }+D_{ij}(P_{bone}^{\max }-P_{nor}^{\max })]\end{array},\forall ij\in S_{ori}$

(16) $(H_{ij}-1)M\le {\widehat{p}}_{ij}^{\left(0\right)}-({\theta }_i^{\left(0\right)}-{\theta }_j^{\left(0\right)})/x_{ij}\le (H_{ij}-1)M,\hspace{1em}\forall ij\in S_{plan}$

(17) $\{\begin{array}{c}{\widehat{p}}_{ij}^{\left(0\right)}\ge (1-{\lambda }_l)[-H_{ij}{P}_{nor}^{\max }-D_{ij}(P_{bone}^{\max }-P_{nor}^{\max })]\\ {\widehat{p}}_{ij}^{\left(0\right)}\le (1-{\lambda }_l)[H_{ij}{P}_{nor}^{\max }+D_{ij}(P_{bone}^{\max }-P_{nor}^{\max })]\end{array},\forall ij\in S_{plan}$

(18) ${\theta }_i^{\min }\le {\theta }_i^{\left(0\right)}\le {\theta }_i^{\max },\hspace{1em}\forall i\in S_N$

(19) $0\le P_{Gi}^{(0)}\le (1-{\lambda }_g)P_{Gi}^{\max },\hspace{1em}\forall i\in S_g$

(20) $0\le \Delta p_i^W\le {\overline{P}}_i^W,\hspace{1em}\forall i\in S_w$

(21) $H_{ij}\ge D_{ij},\hspace{1em}\forall ij\in S_{plan}$

(22) $P_{Gi}^{(d)}-P_{L_i}^{(d)}-{\displaystyle \sum _{ij\in S_{all}}{p}_{ij}^{(d)}}=0,\forall i\in S_N$

(23) $(D_{ij}-1)M_1\le p_{ij}^{\left(d\right)}-({\theta }_i^{\left(d\right)}-{\theta }_j^{\left(d\right)})/x_{ij}\le (1-D_{ij})M_1,\hspace{1em}\forall ij\in S_{all}$

(24) $-D_{ij}{P}_{bone}^{\max }\le p_{ij}^{\left(d\right)}\le D_{ij}{P}_{bone}^{\max },\hspace{1em}\forall ij\in S_{all}$

(25) ${\theta }_i^{\min }\le {\theta }_i^{\left(d\right)}\le {\theta }_i^{\max },\hspace{1em}\forall i\in S_N$

(26) $0\le P_{Gi}^{(d)}\le P_{Gi}^{\max },\hspace{1em}\forall i\in S_g$

(27) $P_{L_i}^{(d)\min }\le P_{L_i}^{(d)}\le P_{L_i}^{(0)},\hspace{1em}\forall i\in S_N$

Equations (13)–(21) are the constraints under normal operation, where $P_{Gi}^{(0)}$ represents the unit output connected to bus i, ${\overline{P}}_i^W$ represents the predicted value of wind power connected to bus i, and $\Delta p_i^W$ represents its corresponding predicted wind power’s abandonment. $P_{L_i}^{(0)}$ represents the load of bus i; $p_{ij}^{(0)}$, ${\widehat{p}}_{ij}^{(0)}$ represent the power flow of the original line and the flow assumption of the optional construction line, respectively; ${\theta }_i^{\left(0\right)}$ and ${\theta }_j^{\left(0\right)}$ represent the voltage phase angle of bus i and bus j, respectively; x_ij represents the reactance of line ij. $P_{Gi}^{\max }$ represents the maximum output of the unit connected to bus i, and λ_g represents its reserve margin. $P_{nor}^{\max }$, $P_{bone}^{\max }$ represent the maximum power flow of the original line and the enhanced line, respectively, and λ_l represents the margin for subsequent regulation. ${\theta }_i^{\max }$, ${\theta }_i^{\min }$ represent the upper and lower limits of the voltage amplitude, respectively; S_g represents the set of buses connected to the generator set; S_w represents the set of buses connected to the wind power; M₁ represents an extremely big number.

Equations (22)–(27) are the constraints under disasters, where $P_{Gi}^{(d)}$ represents the guaranteed output of unit at bus i; $P_{L_i}^{(d)}$ represents the guaranteed load at bus i; and ${\theta }_i^{\left(d\right)}$, ${\theta }_j^{\left(d\right)}$ are the voltage phase angle values of bus i and bus j, respectively; $P_{L_i}^{(d)\min }$ represents the minimum load requirement necessary to supply at bus i, which represents the critical load; S_all represents the set of all lines in the system.

Equation (13) is the power balance constraint under normal operation. Equation (14) is the power flow balance equation. Equation (15) is the capacity constraint of the original line. The expression of the upper and lower bounds takes into account the fact that the capacity will be increased after enhancement. Equations (16) and (17) use the “big-M” method [18] to separate the 0–1 decision variables from the continuous variables in the flow equation of the optional construction lines. The nonlinear limitation of the original constraint is effectively solved. Here, ${\widehat{p}}_{ij}^{\left(0\right)}$ for the optional construction lines is introduced. When H_ij = 1, which means a new line ij is selected, Equation (16) is converted into a flow equation; when H_ij = 0 and D_ij = 0, which means line ij is not built and known from Equation (17), ${\widehat{p}}_{ij}^{\left(0\right)}=0$, and Equation (16) is converted into a slack constraint. Equation (18) is the bus voltage phase angle constraint under normal operation. Equation (19) is the generator output constraint with a margin to cope with the uncertainty. Equation (20) is the disposal limit constraint for predicted wind power. Equation (21) is the enhanced constraint for new lines. Only optional lines that have been newly constructed can be further enhanced for construction.

Equation (22) is the power balance constraint under disasters. According to the actual situation, it can be generally considered that all renewable energy sources such as wind power are shut down under extreme disasters. Considering the more extreme cases, all unenhanced lines are withdrawn from operation. Equations (23) and (24) are the power flow equations of the enhanced lines by using the “big-M” method. Equation (25) is the bus voltage phase angle constraint under disasters. Equation (26) is the output constraint of generators under disasters. Equation (27) is the load guarantee constraint.

In addition, the skeleton network should also satisfy the connectivity constraint, which means the critical loads can be connected through the skeleton network without uncoupling. Here, the “virtual flow” method [29] is used to judge the connectivity: for example, consider a skeleton network given in Figure 3a, and construct a “virtual power supply” g external to a bus a. If a bus in the skeleton network contains a critical load, it is assumed that it has 1 unit of “virtual load” (indicated by the red arrow) and that the “virtual load” can only be supplied by g. Since there is only one “virtual power supply”, it is analogous to the real power flow. Each “virtual load” and g must be connected to each other, and there can be no islands (if not connected, as in Figure 3b, bus i, k, and n will be no supply).

Each branch’s “virtual flow” exists independently of the real power flow. Therefore, if the planned skeleton network is connected, the constraints should be satisfied for each bus i that contains critical loads:

(28)${\displaystyle \sum _{ij\in S_{all}}{f}_{ij}=1}$

(29)$-D_{ij}{f}_{ij}^{max}\le f_{ij}\le D_{ij}{f}_{ij}^{max},ij\in S_{all}$

where f_ij represents the “virtual flow” flowing through the branch ij in the skeleton network; $f_{ij}^{max}$ is the upper limit of the “virtual flow”, which can be generally set as the number of total buses. Equation (28) is the flow balance constraint. Equation (29) is the upper limit constraint of “virtual flow”, and when D_ij = 0, f_ij = 0, the ordinary lines can be distinguished from the enhanced lines.

4.2.2. Lower-Level Model

The lower planning model is an operational calibration level to check the ability to withstand uncertainties for the planning network scheme, which is passed down from the upper model. And the uncertainties considered are mainly caused by the volatility of wind power. The wind power can be fitted according to the GMM described in Section 2. The lower-level planning model further considers whether the network can carry the changes caused by the uncertainty, after receiving the grid planning scheme and relevant operating parameters from the upper-level model. If necessary, it is required to increase the amount of abandoned wind power. Therefore, the objective function of the lower model is to increase the minimum amount of abandoned wind power:

(30)$\min f(\Delta p_{ad}^W)={\displaystyle \sum _{j\in S_w}\Delta p_{j.ad}^W}$

where $\Delta p_{j.ad}^W$ represents the additional abandoned wind power at bus j, which is based on the abandoned wind power determined at the upper level.

For the aim to cope with the uncertainty, considering the common imitation control strategy for units in engineering, the actual output of the generators set depending on the actual situation can be expressed as:

(31)$P_{gi}=P_{Gi.set}-{\alpha }_i{\displaystyle \sum _{j\in S_w}[{\tilde{P}}_j^W-(\Delta p_{j.set}^W+\Delta p_{j.ad}^W)-({\overline{P}}_j^W-\Delta p_{j.set}^W)]}=P_{Gi.set}-{\alpha }_i{\displaystyle \sum _{j\in S_w}({\tilde{\epsilon }}_j^W-\Delta p_{j.ad}^W)},\forall i\in S_g$

(32)${\displaystyle \sum _{i\in S_g}{\alpha }_i=1}$

where P_Gi.set represents the predicted power output of the unit obtained from the upper-level model; $\Delta p_{j.set}^W$ represents the predicted wind power’s abandonment obtained from the upper-level optimization; α_i represents the power allocation factor of the i_th unit.

For the network structure given in the upper-level, the branch ij’s power flow p_ij can be obtained as:

(33)$p_{ij}={\displaystyle \sum _{o\in S_g}{G}_{go.ij}{P}_{go}}+{\displaystyle \sum _{w\in S_w}{G}_{w.ij}({\tilde{P}}_w^W-\Delta p_{w.set}^W-\Delta p_{w.ad}^W)}+{\displaystyle \sum _{d\in S_N}{G}_{d.ij}{P}_{Ld}^{(0)}}$

where G_go.ij, G_w.ij, and G_d.ij are the active power transfer distribution factors of the branch ij to the units, wind power, and loads, respectively.

The operational constraints that should be satisfied by the lower-level model are

s.t.

(34)${\displaystyle \sum _{i\in S_g}{P}_{gi}}+{\displaystyle \sum _{j\in S_w}({\tilde{P}}_j^W-\Delta p_{j.set}^W-\Delta p_{j.ad}^W)}={\displaystyle \sum _{i\in S_N}{P}_{Li}^{(0)}}$

(35)$-P_{nor}^{\max }-D_{ij}^{set}(P_{bone}^{\max }-P_{nor}^{\max })\le p_{ij}\le P_{nor}^{\max }+D_{ij}^{set}(P_{bone}^{\max }-P_{nor}^{\max }),\forall ij\in S_{ori}$

(36)$-H_{ij}^{set}{P}_{nor}^{\max }-D_{ij}^{set}(P_{bone}^{\max }-P_{nor}^{\max })\le p_{ij}\le H_{ij}^{set}{P}_{nor}^{\max }+D_{ij}^{set}(P_{bone}^{\max }-P_{nor}^{\max }),\forall ij\in S_{plan}$

(37)$0\le P_{gi}\le P_{Gi}^{\max },\hspace{1em}\forall i\in S_g$

where $H_{ij}^{set}$ and $D_{ij}^{set}$ are the planning results of the new ordinary line and the new enhanced construction line obtained from the upper level, respectively.

Equation (34) is the power balance constraint of the system. Equation (35) is the capacity constraint of the original lines. Equation (36) is the capacity constraint of the optional lines. Equation (37) is the power output constraint of the generators.

4.2.3. The Relationships between Upper and Lower-Level Model

The network planning scheme $H_{ij}^{set}$, $D_{ij}^{set}$ and the corresponding operating parameters P_Gi.set, $\Delta p_{j.set}^W$ are obtained from the upper-level optimization. Based on this, the lower-level model further finds the changes triggered by the uncertainty to carry wind power. And the additional amount of abandoned wind power $\Delta p_{ad}^W$ can be obtained and feedback to the upper level to improve the objective function. The total objective function value of this network structure scheme is obtained. Finally, the optimal network structure solution can be selected that considers both wind power uncertainty and disasters.

5. Convert and Solve Chance Constrained Model Based on Convex Relaxation

5.1. Chance Constraint

The main manifestation of wind power’s uncertainty is that the output is in the form of a probability distribution, which is a set of random variables. And the chance constraint can better characterize the uncertainty brought by random variables [25], which can be generally expressed as:

(38)$\mathrm{Pr}\{y(x,\tilde{\epsilon })\ge 0\}\ge 1-\eta$

where x is the decision variable; $\tilde{\epsilon }$ is the random variable; and η is the violation probability of this chance constraint. Considering the influence of uncertainty factors in the power system, it is generally preferred to control the operating constraints within the controllable range. When η = 0, the confidence level is 1, indicating that the operational constraints must be strictly satisfied for all possible situations. Although it can greatly improve the security, it will inevitably also cause the waste of resources.

Therefore, considering the impact of wind power’s uncertainty, Equations (35) and (36) are written uniformly in the form of the following chance constraint:

(39)$\mathrm{Pr}\{p_{ij}+P_{ij}^{\max }\ge 0\}\ge 1-\eta ,\forall ij\in S_{all}$

(40)$\mathrm{Pr}\{-p_{ij}+P_{ij}^{\max }\ge 0\}\ge 1-\eta ,\forall ij\in S_{all}$

where $P_{ij}^{\max }$ is the maximum tidal limit of the known line ij.

5.2. The Reduction of Model Size

5.2.1. The Influence Area of Wind Power

Given that the maximum power flow margin is left in the upper model, it is not necessary to impose the mandatory restriction on the power flow limit constraint of each line. Therefore, the concept of the influence area of wind farm is defined, aiming to select the set of lines that are affected by the fluctuation of wind power. The reference coefficient is introduced as shown in Equation (41):

(41)${\beta }_{w.ij}=\left|\frac{p_{ij}^{w\max }-p_{ij}^0}{p_{ij}^0}\right|\times 100\%$

where $p_{ij}^0$ represents the power flow of line ij when all wind power is 0; $p_{ij}^{w\max }$ represents the power flow when the wind farm bus w is of full power and the rest of the wind farms are out of power. And when β_w.ij exceeds a certain threshold τ, the power flow of line ij can be considered more susceptible to the bus w, and it will be included in the influence area of bus w. Eventually, the total influence area of all wind farms can be obtained.

5.2.2. WARD Equivalent

The lines that are more influenced by wind power tend to be concentrated in partial lines. Therefore, the above discussion narrows the range of lines considered by the constraint. Furthermore, in dealing with the constraints of Equations (39) and (40), the WARD equivalent method [30] is used to equate the system structure and reduce the size.

Figure 4 shows a two-connection region system before and after WARD equivalence. The extranet is the region to be equated, and the intranet is the analytical region. The two are connected by boundary nodes Y_i…Y_j and outer network boundary nodes Y_m…Y_n. The WARD equivalent is achieved by eliminating the external nodes and performing the equivalent process: by injecting currents at the boundary nodes and by simulating the action of the extranet on the intranet through the equivalence-to-ground branches of the boundary nodes and the equivalence branches between the boundary nodes.

The DC equation of the system before WARD equivalence can be expressed as:

(42)$\left[\begin{array}{ll}{\mathit{B}}_{\mathbf{EE}}^{(0)}\hfill & {\mathit{B}}_{\mathbf{EI}}^{(0)}\hfill \\ {\mathit{B}}_{\mathbf{IE}}^{(0)}\hfill & {\mathit{B}}_{\mathbf{II}}^{(0)}\hfill \end{array}\right]\left[\begin{array}{l}{\mathit{\theta }}_{\mathbf{E}}\hfill \\ {\mathit{\theta }}_{\mathbf{I}}\hfill \end{array}\right]=\left[\begin{array}{l}{\mathit{P}}_{\mathbf{E}}^{(0)}\hfill \\ {\mathit{P}}_{\mathbf{I}}^{(0)}\hfill \end{array}\right]$

where ${\mathit{B}}_{\mathbf{II}}^{(0)}$ represents the initial chunked susceptance matrix of the intranet; ${\mathit{B}}_{\mathbf{EE}}^{(0)}$ represents the initial chunked susceptance matrix of the extranet nodes; ${\mathit{B}}_{\mathbf{EI}}^{(0)}$ represents the initial chunked susceptance matrix of the extranet to the intranet; ${\mathit{B}}_{\mathbf{IE}}^{(0)}$ represents the initial chunked susceptance matrix of the intranet to the extranet.

After eliminating the external nodes, the power equation is:

(43)${\mathit{B}}_{\mathbf{I}}{}^{\prime }{\mathit{\theta }}_{\mathit{I}}{}^{\prime }={\mathit{P}}_{\mathit{I}}{}^{\prime }$

where ${\mathit{B}}_{\mathit{I}}={\mathit{P}}_{\mathbf{I}}^{(0)}-{\mathit{B}}_{\mathbf{IE}}^{(0)}{\mathit{B}}_{\mathbf{EE}}^{(0)-1}{\mathit{B}}_{\mathbf{EI}}^{(0)}$; ${\mathit{P}}_{\mathit{I}}={\mathit{P}}_{\mathrm{I}}^{(0)}-{\mathit{B}}_{\mathbf{IE}}^{(0)}{\mathit{B}}_{\mathbf{EE}}^{(0)-1}{\mathit{P}}_{\mathrm{E}}^{(0)}$.The key to the WARD equivalent is the determination of the equivalent range. Here, the aim is to select the potential lines that are susceptible to wind power. Therefore, the internal nodes of WARD equivalent include:

• All nodes connected with wind farms;

• All nodes at both ends of the lines within the wind farms’ influence area;

• The nodes at both ends of the optional construction lines.

5.3. Chance Constrained Transformation Based on Convex Relaxation

5.3.1. Convex Relaxation Method

The convex relaxation method converts chance constrained inequalities containing random variables into probability mean inequalities [31]. Thus, the uncertainty of the random variable is characterized in the form of a probability inequality.

For the internal inequality in Equation (38), without loss of universality, it is written in the form of a radiative combination of random variables:

(44)$y(x,\tilde{\epsilon })=y_0(x)+{\displaystyle \sum _{i=1}^N{a}_i{\tilde{\epsilon }}_i\ge 0}$

where a_i is the coefficient before the random variable; y₀(x) is the rest of the equation.

Both decision variables and random variables have upper and lower bounds, so $y(x,\tilde{\epsilon })$ must have upper and lower bounds. It may be useful to set its lower bound as −L (L > 0), so that $g(x,\tilde{\epsilon })=y(x,\tilde{\epsilon })/L+1$, and introduce a certain monotonic non-negative concave function φ; then, Equation (38) can be transformed into

(45)$\mathrm{Pr}\{\phi (g(x,\tilde{\epsilon }))\ge \phi (1)\}\ge 1-\eta$

Using Markov’s inequality [32], it is further obtained that:

(46)$E[\phi (g(x,\tilde{\epsilon }))]\ge (1-\eta )\phi (1)$

where E[·] denotes taking the mean value. Regarding the selection of functions, the literature gives a class of functions with good relaxation effects:

(47)$\phi (z)=1-e^{-z}$

Jointly union Equations (44), (46), and (47), the following can be obtained:

(48)${\displaystyle \sum _{i=1}^N\ln E(e^{-a_i{\tilde{\epsilon }}_i/L})}-\frac{y_0(x)}{L}\le \ln (1-\eta +e\eta )$

If ${\tilde{\epsilon }}_i$ is in the form of the GMM as in Equation (1), it can be deduced (see the Appendix A for the procedure) that:

(49)$E(e^{-a_i{\tilde{\epsilon }}_i/L})={\displaystyle \sum _{m=1}^M{\omega }_{m,i}{e}^{\frac{a_i^2{\sigma }_{m,i}^2}{2L^2}-\frac{a_i{\mu }_{m,i}}{L}}}$

Therefore, Equation (38) can be relaxed into the inequality:

(50)${\displaystyle \sum _{i=1}^N\ln {\displaystyle \sum _{m=1}^M{\omega }_{m,i}{e}^{\frac{a_i^2{\sigma }_{m,i}^2}{2L^2}-\frac{a_i{\mu }_{m,i}}{L}}}}-\frac{y_0(x)}{L}\le \ln (1-\eta +e\eta )$

Regarding the selection of the lower bound L, it is necessary to keep iterating to improve the relaxation effect by increasing the accuracy of the taken values. The iterative process is discussed below:

Assuming a_i of Equation (44) is greater than 0, if not, the probability distribution of ${\tilde{\epsilon }}_i$ can be made equivalent by taking negative. According to Equation (44), it may be assumed to have been written in the form of GMM. Each ${\tilde{\epsilon }}_i$ satisfies a single Gaussian distribution, which is not contradictory to the original problem. Considering the actual situation, ${\tilde{\epsilon }}_i$ can be approximated the lower bound of μ_m,I − 3σ_m,i. And assume that the initial lower bound for y₀(x) is l₀. Therefore, the initial L₀ of L’s lower bound can be made as follows:

(51)$L_0=-[l_0+{\displaystyle \sum _{i=1}^N{a}_i({\mu }_{m,i}-}3{\sigma }_{m,i})]$

From the relaxation of Equation (50), the following is obtained for the k_th iteration:

(52)$y_k(x)\ge {\displaystyle \sum _{i=1}^N(\frac{a_i^2{\sigma }_{m,i}^2}{2L_k}-a_i{\mu }_{m,i})}-L_k\ln (1-\eta +e\eta )$

Define s = ln(1 − η + eη); then, l₀ of the k_th iteration can be taken as:

(53)$l_k={\displaystyle \sum _{i=1}^N(\frac{a_i^2{\sigma }_{m,i}^2}{2L_k}-a_i{\mu }_{m,i})}-L_{k}s$

According to Equation (51), L_k for the k_th iteration can be deduced as:

(54)$L_k=-[{\displaystyle \sum _{i=1}^N(\frac{a_i^2{\sigma }_{m,i}^2}{2L_{k-1}}-3a_i{\sigma }_{m,i})}-L_{k-1}s]$

It is known that Equation (54) is a recursive formula for L, which can be made to converge to the steady state point L_x after constant recursion, satisfying:

(55)$(1-s)L_x^2-L_x{\displaystyle \sum _{i=1}^{N}3a_i{\sigma }_{m,i}}+{\displaystyle \sum _{i=1}^N{a}_i^2{\sigma }_{m,i}^2}=0$

Convergence is the solution of Equation (55). To ensure that L_x > 0, take the larger solution, that is, the sought L:

(56)$L=\frac{{\displaystyle \sum _{i=1}^N\left(3a_i{\sigma }_{m,i}\right)}+\sqrt{{\left({\displaystyle \sum _{i=1}^{N}3}{a}_i{\sigma }_{m,i}\right)}^2-2(1-s){\displaystyle \sum _{i=1}^N{a}_i^2{\sigma }_{m,i}^2}}}{2(1-s)}$

5.3.2. Transformation of Chance Constrained Planning Model

The p_ij in Equations (39) and (40) is affected by the uncertainty of wind power. To quantify this uncertainty, p_ij needs to be expressed in terms of the related random variables.

After determining the internal node objects to be WARD equivalent in Section 5.2, from Equation (36), it can be obtained:

(57)${\mathit{\theta }}_{\mathit{I}}{}^{\prime }={({\mathit{B}}_{\mathit{I}}{}^{\prime })}^{-1}{\mathit{P}}_{\mathit{I}}{}^{\prime }={\mathit{G}}_{\mathit{I}}{}^{\prime }{\mathit{P}}_{\mathit{I}}{}^{\prime }$

where ${\mathit{G}}_{\mathit{I}}{}^{\prime }={({\mathit{B}}_{\mathit{I}}{}^{\prime })}^{-1}$, combining Equations (31)–(33), it can be obtained:

(58)$\begin{array}{ll}{p}_{ij}=\frac{1}{x_{ij}}({\theta }_i-{\theta }_j)=\frac{1}{x_{ij}}{\displaystyle \sum _{k=1}^n(G_{ik}-G_{jk})P_k}& =\frac{1}{x_{ij}}{\displaystyle \sum _{k=1}^n(G_{ik}-G_{jk})\{[P_{Gk.set}-{\alpha }_k{\displaystyle \sum _{w\in S_w}({\tilde{\epsilon }}_w^W-\Delta p_{w.ad}^W)}]-P_{L_i}^{(0)}\}}\\ & +\frac{1}{x_{ij}}{\displaystyle \sum _{w=1}^W(G_{iw}-G_{jw})({\overline{P}}_w^W+{\tilde{\epsilon }}_w^W-\Delta p_{w.set}^W-\Delta p_{w.ad}^W)}\end{array}$

where G_ik and G_jk are the corresponding elements of the n-dimensional matrix. To facilitate the subsequent calculation of the expected value, note that C_ij will be the constant part of Equation (58), S_ij.w is the coefficient before the variable to be optimized in Equation (58), and R_ij.w is the coefficient before the random variable in Equation (58), then:

(59)$C_{ij}=\frac{1}{x_{ij}}{\displaystyle \sum _{k=1}^n(G_{ik}-G_{jk})(P_{Gk.set}-P_{L_i}^{(0)})}+\frac{1}{x_{ij}}{\displaystyle \sum _{w=1}^W(G_{iw}-G_{jw})({\overline{P}}_w^W-\Delta p_{w.set}^W)}$

(60)$S_{ij.w}=\frac{1}{x_{ij}}{\displaystyle \sum _{k=1}^n(G_{ik}-G_{jk}){\alpha }_k}+\frac{1}{x_{ij}}{\displaystyle \sum _{w=1}^W(-G_{iw}+G_{jw})}$

(61)$R_{ij.w}=\frac{1}{x_{ij}}{\displaystyle \sum _{k=1}^n(-G_{ik}+G_{jk})}+\frac{1}{x_{ij}}{\displaystyle \sum _{w=1}^W(G_{iw}-G_{jw})}$

Then, it is obtained that:

(62)$p_{ij}=C_{ij}+{\displaystyle \sum _{w\in S_w}{S}_{ij.w}\Delta p_{w.ad}^W}+{\displaystyle \sum _{w\in S_w}{R}_{ij.w}{\tilde{\epsilon }}_w^W}$

Further, the chance constraint of Equation (39) can be transformed based on the above convex relaxation method: note

(63)$y_1(\Delta p_{ad}^W,{\tilde{\epsilon }}^W)=p_{ij}+P_{ij}^{\max }=P_{ij}^{\max }+C_{ij}+{\displaystyle \sum _{w\in S_w}{S}_{ij.w}\Delta p_{w.ad}^W}+{\displaystyle \sum _{w\in S_w}{R}_{ij.w}{\tilde{\epsilon }}_w^W}$

And assume that the lower bound of $y_1(\Delta p_{ad}^W,{\tilde{\epsilon }}^W)$ is L₁; then Equation (39) is converted to:

(64)${\displaystyle \sum _{w\in S_w}\ln {\displaystyle \sum _{m=1}^M{\omega }_{m,w}{e}^{\frac{R_{ij.w}^2{\sigma }_{m,w}^2}{2L_1^2}-\frac{R_{ij.w}{\mu }_{m,w}}{L_1}}}}-\frac{1}{L_1}(P_{ij}^{\max }+C_{ij}+{\displaystyle \sum _{w\in S_w}{S}_{ij.w}\Delta p_{w.ad}^W})\le s$

where:

(65)$L_1=\frac{{\displaystyle \sum _{w\in S_w}(3R_{ij.w}{\sigma }_{m,w})}+\sqrt{{\left({\displaystyle \sum _{w\in S_w}3R_{ij.w}{\sigma }_{m,w}}\right)}^2-2(1-s){\displaystyle \sum _{w\in S_w}{R}_{ij.w}^2{}_{ij.w}{\sigma }_{m,w}^2}}}{2(1-s)}$

Similarly, the chance constraint of Equation (40) can be transformed:

(66)$y_2(\Delta p_{ad}^W,{\tilde{\epsilon }}^W)=-p_{ij}+P_{ij}^{\max }=P_{ij}^{\max }-C_{ij}-{\displaystyle \sum _{w\in S_w}{S}_{ij.w}\Delta p_{w.ad}^W}-{\displaystyle \sum _{w\in S_w}{R}_{ij.w}{\tilde{\epsilon }}_w^W}$

Assume that the lower bound of $y_2(\Delta p_{ad}^W,{\tilde{\epsilon }}^W)$ is L₂; then Equation (40) is converted to:

(67)${\displaystyle \sum _{w\in S_w}\ln {\displaystyle \sum _{m=1}^M{\omega }_{m,w}{e}^{\frac{R_{ij.w}^2{\sigma }_{m,w}^2}{2L_2^2}+\frac{R_{ij.w}(-{\mu }_{m,w})}{L_2}}}}-\frac{1}{L_2}(P_{ij}^{\max }-C_{ij}-{\displaystyle \sum _{w\in S_w}{S}_{ij.w}\Delta p_{w.ad}^W})\le s$

where:

(68)$L_2=\frac{{\displaystyle \sum _{w\in S_w}(3R_{ij.w}{\sigma }_{m,w})}+\sqrt{{\left({\displaystyle \sum _{w\in S_w}3R_{ij.w}{\sigma }_{m,w}}\right)}^2-2(1-s){\displaystyle \sum _{w\in S_w}{R}_{ij.w}^2{}_{ij.w}{\sigma }_{m,w}^2}}}{2(1-s)}$

In summary, the integrated the network of transmission planning model can then be restated as follows:

Upper-level model:

(69)$\min F={\displaystyle \sum _{t\in S_{plan}}[C_{nor}{H}_t{L}_t}+C_{bone}{D}_t{L}_t]+{\displaystyle \sum _{t\in S_{ori}}{C}_{bone}{D}_t{L}_t}+{\displaystyle \sum _{j\in {\varphi }_w}\alpha \Delta p_j^W}-{\chi }_L{\displaystyle \sum _{i\in S_N}{P}_{L_i}^{(d)}}+\delta f(\Delta p_{ad}^W)$

s.t Equations (13)–(29)

Lower-level model:

(70)$\min f(\Delta p_{ad}^W)={\displaystyle \sum _{j\in S_w}\Delta p_{j.ad}^W}$

s.t Equations (31), (32), (34), (37), (59)–(61), (64), (65), (67) and (68).

Figure 5 gives the flow chart for solving the integrated network planning problem.

6. Case Study

6.1. Parameters and Operating Environment Configuration

For the aim to verify the effectiveness of the proposed integrated transmission network planning method, simulation tests are conducted in the IEEE118 bus system. The topology of the IEEE118 bus system is shown in Figure 6, which contains 118 buses and 186 lines (represented by black thin solid lines) and 15 lines that can be selected for construction (represented by blue dashed lines, whose resistance and reactance information is shown in Appendix B). There are 54 buses connected with generators, represented by triangles. And the rest are substation buses, represented by circles. The white triangles and circles indicate these buses aren’t connected with loads. Forty-nine buses are marked in red to indicate that they contain critical loads. Detailed critical load information is given in Appendix B. Here it is assumed that bus 6, 56, and 85 are connected to wind farms.

The length of line is calculated by converting the reactance to 0.320 $\mathrm{\Omega }$/km, and the investment cost C_nor and the cost C_bone are taken as 24,900 dollar/km and 25,000 dollar/km, respectively. The penalty factor for wind power’s abandonment α and the penalty factor δ are both taken as 10⁴. The maximum allowable power flow for ordinary lines $P_{nor}^{\max }$ and enhanced lines $P_{bone}^{\max }$ are taken as 100 MW and 200 MW, respectively. In addition, the application of transmission cables is becoming more and more widespread, and considering that cables are usually buried underground or in tunnels, they have higher disaster resistance compared to overhead lines. And they are generally regarded as the base lines of the skeleton network. The specific information of cable lines is shown in Appendix B.

The program for solving the bi-level model of the transmission network planning is written by C++, and the commercial solver Cplex is called to perform the solution. When solving the lower-level model, the solution pool parameters of Cplex are set as follows: intensity is set to 4; capacity is set to 200; relgap is set to 0.1 and replaced to 2. When solving the upper-level model, the parameters are used in their default configuration. The simulation results are conducted on a PC with Intel Core i7-8550U CPU @1.99GHZ with 8.00GB RAM.

6.2. IEEE 118 Bus System

6.2.1. Quantitative Analysis of Wind Power’s Uncertainty

The historical output data of three wind farms connected to the system were obtained from the open-source data provided by the elia group [33]. The predicted wind power ${\overline{P}}_6^W$, ${\overline{P}}_{56}^W$, and ${\overline{P}}_{85}^W$ are taken as 400 MW, 400 MW, and 600 MW, respectively. Based on the historical output data and AIC, the number of better results for fitting ${\tilde{\epsilon }}_6^W$, ${\tilde{\epsilon }}_{56}^W$, and ${\tilde{\epsilon }}_{85}^W$ are 6, 7, and 7, respectively. The GMM is used to fit, and Figure 7 gives the fitting results for the deviations.

It can be seen from Figure 7 that each figure contains the frequency distribution and probability density function (PDF) of the predicted deviation of the wind power. And the PDF curve is a good fit for the deviation distribution trend. On the other hand, it can also be seen that the output deviation has a single-peaked nature, and the larger the output deviation is, the lower the probability of occurrence is. This is compatible with the factual situation and indicates that the GMM can better describe the uncertainty of the wind power. GMM will also provide a foundation for subsequent convex relaxation.

6.2.2. The Integrated Network Planning without Considering Wind Power’s Uncertainty

When wind power’s uncertainty is not considered, only the upper-level model of the integrated network planning described in Section 4.2.1 is analyzed. The upper model is solved directly by removing $f(\Delta p_{ad}^W)$ which is returned from the lower level. The schematic diagram of the integrated network planning without considering wind power’s uncertainty is shown in Figure 8. The red thick solid line indicates the enhanced construction of original lines, and the blue thick solid line indicates the enhanced optional construction lines, and the two together constitute the skeleton network of the system.

Table 1 gives the specific network planning scheme of the transmission. If the fluctuation of wind power is not considered, combined with Figure 8, it can be seen that 68 lines are enhanced and none of lines are selected as ordinary lines in the original lines. Moreover, four lines are chosen to be built as enhanced lines, and none of them are as ordinary lines in the optional lines. It means that the increase of original lines and capacity of enhanced lines is able to meet the requirements of wind power access. All ordinary lines will be withdrawn from operation due to the disasters, where more consideration is given to the critical loads. For example, bus 16 contains some critical loads, and it has two options to connect to the skeleton network: enhance original line 16–17 or build a new enhanced line 6–16. However, the length of line 16–17 is more than seven times that of line 6–16, and the new construction of the latter is more economical in view of investment cost.

It can be observed that bus 77 and bus 100 contain the larger critical loads and they are with a high critical level. Therefore, a larger number of enhanced lines is needed to be connected in order to secure the critical load supply. It’s noted that lines 30–38, 37–38, 65–68, and 38–68 have no significant loads at either end of the lines but these lines are also selected for enhanced construction as the key connecting hubs of the network frame. In summary, it shows that the results of the planned network are in line with the actual requirements.

All 72 enhanced lines together form a topologically connected skeleton network, indicating that the implementation of “virtual flows” is effective. Moreover, all critical loads are covered in the skeleton network, which means that only these 72 enhanced lines are needed to ensure the supply of critical loads under disasters. On the one hand, no loops are formed in the skeleton network, and if loops appear, it must mean that the existence of a line in the loop is a non-essential enhancement. On the other hand, the engineering amount of line differential enhancement is only 36.32% compared to the full line enhancement, both of which indicate the better economy of the planned network frame.

6.2.3. The Integrated Network Planning Considering Wind Power’s Uncertainty

Here, the threshold value of wind farms’ influence area τ is taken as 20%. After the calculation, then there are 23 lines in the influence area of the wind farms, and there are 31 internal buses after WARD equivalent, accounting for 26.27% of all buses. Subsequently, only the lines between these buses need to be examined, which can effectively reduce the scale of the calculation. When the probability of chance constraint violation η is taken as 5%, the result of the proposed integrated network planning model considering wind power uncertainty is shown in Figure 9.

In Figure 9 with the red, thick solid line for enhanced construction of original lines, the blue thick solid line indicates the new enhanced lines; the blue thin solid line indicates the new ordinary lines; the specific planning scheme is shown in Table 1. Sixty-six lines are selected to be enhanced in original lines. And four lines are chosen to be built as enhanced lines, and three are chosen to be built as ordinary lines among the optional construction lines. The green dashed circles indicate the changed construction range compared with Figure 8. It can be seen that the changed parts are mainly concentrated near the wind farm buses, which indicates that the area around the buses connected with wind farms is more affected after considering the wind power’s uncertainty, while the ones farther away from these buses are less affected. On the other hand, the new enhanced lines 6–11, 6–16, 77–83, and 85–90 and the new ordinary lines 52–54, 53–56, and 83–88 are located around the wind farms. This shows that new lines are needed to cope with the fluctuating influence after considering the fluctuating influence of wind power and also makes a small correction to the skeleton network under disasters. It effectively proves that the network planning scheme has practical significance.

The robustness of the planned network scheme is further analyzed: according to the previous GMM fitting results, 1000 scenarios are randomly sampled for all wind farms. Then, they are tested for whether the system satisfies the operational constraints by substituting into the planned network structure. Note that N_pass is the number of scenarios that satisfy the running constraint, and the pass ratio P_Ra is:

(71)$P_{Ra}=N_{pass}/1000$

The scenarios’ pass ratio under different chance constraint violation probabilities of η is given in Figure 10, where the horizontal coordinate is the threshold value of wind farms’ influence area τ and the vertical coordinate is P_Ra. From the graph, it can be seen that P_Ra decreases as the value of τ increases. This is due to the fact that when the threshold τ is larger, the range of wind farms’ influence area is smaller, leading to a smaller scope of investigation, which may ignore some lines where their power flow cross limits. It can be seen from the figure that when τ is taken as 20%, P_Ra is already better for different chance constraint violation probabilities η. On the other hand, when η is larger, P_Ra is smaller. This is because the larger the value of η, the lower the confidence probability that the chance constraint holds, and the greater the number of scenarios that the planned grid cannot withstand wind power’s uncertainty. However, it can be seen from the figure that when τ is less than 20%, P_Ra is larger than the set confidence level 1 − η, which indicates that the threshold value τ taken in this paper is reasonable. Thus, the effectiveness of the planned grid structure is proved.

7. Conclusions

An integrated network planning method for transmission is proposed in this paper, which comprehensively considers wind power’s uncertainty and load security under disasters. The method facilitates the solution of the model by building up a bi-level model with chance constraints, and the uncertainty of wind power is accurately expressed by using GMM. And then a convex relaxation method is constructed for chance constraints suitable for GMM. Furthermore, the chance constraints of this model are converted into probabilistic inequalities by using this convex relaxation method. The simulation results obtained in the IEEE118 bus system and subsequent scenario tests demonstrate the effectiveness of the proposed integrated network planning model. It is shown that the network planning scheme can withstand wind power’s uncertainty and ensure the supply of critical loads under disasters.

In future work, the correlation between the output for different wind farms can be further explored, and how this correlation influences the integrated network planning model for the transmission grid can be discussed.

Author Contributions

Conceptualization, Y.S. and R.G.; methodology, R.G. and Y.S.; software, R.G. and Y.S; validation, Y.S. and Z.Y.; formal analysis, Y.T. and Y.L.; investigation, Y.S. and Z.Y.; writing—original draft preparation, Y.S., Z.Y., Y.T. and Y.L.; writing—review and editing, Y.S. and Z.Y.; supervision, R.G., Y.T. and Y.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Publicly available datasets were analyzed in Section 6 of this study. This data can be found here: https://matpower.org/docs/ref/matpower7.1/menu7.1.html (accessed on 25 June 2023). https://www.elia.be/en/grid-data/power-generation/wind-power-generation (accessed on 28 June 2023).

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Diagram of skeleton network against disasters and the critical loads.

Figure 2. The framework of integrated network planning for the transmission.

Figure 3. Schematic diagram of “virtual flow” method. (a) contain connectivity (b) not contain connectivity.

Figure 4. A two-connection region system before and after WARD equivalent.

Figure 5. The flow chart for solving the integrated network planning problem for transmission.

Figure 6. The topology diagram of IEEE118 bus system.

Figure 7. The deviation fitting for wind power prediction based on GMM. (a) Deviation fitting of [Forumla omitted. See PDF.] (b) Deviation fitting of [Forumla omitted. See PDF.]. (c) Deviation fitting of [Forumla omitted. See PDF.].

Figure 8. The integrated network planning without considering wind power’s uncertainty.

Figure 9. The integrated network planning considering wind power’s uncertainty.

Figure 10. The scenarios’ pass ratio under different chance constraint violation probability η.

Table 1

Specific network planning scheme of the transmission under two situations.

Situation	Construction Scope	Title 3
Without considering wind power’s uncertainty	original lines	Enhance:1–3,3–5,4–5,4–11,5–6,6–7,5–8,8–30,15–19,15–33,15–17,17–30,17–31,29–31,31–32,32–114,114–115,27–115,30–38,37–38,34–36,34–37,37–39,39–40,40–42,54–55,54–56,55–59,59–63,63–64,60–61,61–62,61–64,64–65,65–66,45–49,46–47,47–49,49–66,65–68,68–81,70–74,74–75,75–118,76–118,76–77,77–78,77–80,80–81,80–98,82–83,83–85,85–88,82–96,95–96,94–95,93–94,92–93,94–100,100–101,100–103,100–104,104–105,105–107,105–108,108–109,109–110
	optional construction lines	Build as ordinary lines: None
		Build as enhanced lines: 6–16,38–68,77–83,85–90
Considering wind power’s uncertainty	original lines	Enhance: 1–3,3–5,4–5,4–11, 6–7, 5–8,8–30,15–19,15–33,15–17,17–30,17–31,29–31,31–32,32–114,114–115,27–115,30–38,37–38,34–36,34–37,37–39,39–40,40–42,54–56,55–56,55–59,59–63,63–64,60–61,61–62,61–64,64–65,65–66,45–49,46–47,47–49,49–66,65–68,68–81,70–74,74–75,75–118,76–118,76–77,77–78,77–80,80–81,80–98,82–83,83–85,85–88,82–96,95–96,94–95,93–94,92–93,94–100,100–101,100–103,100–104,104–105,105–107,105–108,108–109,109–110
	optional construction lines	Build as ordinary lines: 52–54,53–56,83–88
		Build as enhanced lines: 6–11,6–16,77–83,85–90

Appendix A

The derivation process of Equation (49):

When ${\tilde{\epsilon }}_i$ is in the form of the GMM forms as Equation (A1), the probability density function of each Gaussian component ϕ_m,i is assumed to be:(A1)$${\varphi }_{m,}\left({\tilde{\epsilon }}_i\mid {\mu }_{m,i},{\sigma }_{m,i}^2\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{m,i}}{e}^{-{(\tilde{\epsilon }-{\mu }_{m,i})}^2/2{\sigma }_{m,i}^2}$$

Then, it can be obtained:(A2)$$\begin{array}{ll}E(e^{-a_i{\tilde{\epsilon }}_i/L})& ={\displaystyle {\int }_{-\infty }^{+\infty }{e}^{-a_i{\tilde{\epsilon }}_i/L}{\displaystyle \sum _{m=1}^M\frac{{\omega }_{m,i}}{\sqrt{2\pi }{\sigma }_{m,i}}{e}^{-{({\tilde{\epsilon }}_i-{\mu }_{m,i})}^2/2{\sigma }_{m,i}^2}}d{\tilde{\epsilon }}_i}\\ & ={\displaystyle \sum _{m=1}^M[\frac{{\omega }_{m,i}}{\sqrt{2\pi }{\sigma }_{m,i}}}{\displaystyle {\int }_{-\infty }^{+\infty }{e}^{-\frac{{\tilde{\epsilon }}_i}{2{\sigma }_{m,i}^2}+\frac{{\tilde{\epsilon }}_i{\mu }_{m,i}}{{\sigma }_{m,i}^2}-\frac{a_i}{L}{\tilde{\epsilon }}_i-\frac{{\mu }_{m,i}^2}{2{\sigma }_{m,i}^2}}d{\tilde{\epsilon }}_i]}\\ & ={\displaystyle \sum _{m=1}^M[\frac{{\omega }_{m,i}}{\sqrt{2\pi }{\sigma }_{m,i}}}{\displaystyle {\int }_{-\infty }^{+\infty }{e}^{-\frac{1}{2{\sigma }_{m,i}^2}{({\tilde{\epsilon }}_i-({\mu }_{m,i}-\frac{a_i}{L}{\sigma }_{m,i}^2))}^2+\frac{a_i^2}{2L^2}{\sigma }_{m,i}^2-\frac{a_i{\mu }_{m,i}}{L}}d{\tilde{\epsilon }}_i]}\\ & ={\displaystyle \sum _{m=1}^M[\frac{{\omega }_{m,i}}{\sqrt{2\pi }{\sigma }_{m,i}}{e}^{\frac{a_i^2}{2L^2}{\sigma }_{m,i}^2-\frac{a_i{\mu }_{m,i}}{L}}}{\displaystyle {\int }_{-\infty }^{+\infty }{e}^{-\frac{1}{2{\sigma }_{m,i}^2}{({\tilde{\epsilon }}_i-({\mu }_{m,i}-\frac{a_i}{L}{\sigma }_{m,i}^2))}^2}d({\tilde{\epsilon }}_i-({\mu }_{m,i}-\frac{a_i}{L}{\sigma }_{m,i}^2)]}\\ & ={\displaystyle \sum _{m=1}^M(\frac{{\omega }_{m,i}}{\sqrt{2\pi }{\sigma }_{m,i}}{e}^{\frac{a_i^2}{2L^2}{\sigma }_{m,i}^2-\frac{a_i{\mu }_{m,i}}{L}}}\cdot \sqrt{2\pi }{\sigma }_{m,i})={\displaystyle \sum _{m=1}^M{\omega }_{m,i}{e}^{\frac{a_i^2{\sigma }_{m,i}^2}{2L^2}-\frac{a_i{\mu }_{m,i}}{L}}}\end{array}$$

Appendix B

Table A1

Resistance and reactance information of optional construction lines.

From Bus	To Bus	R (p.u.)	X (p.u.)
Bus-6	Bus-11	0.00862	0.2570
Bus-6	Bus-16	0.00342	0.0514
Bus-7	Bus-11	0.02225	0.0731
Bus-12	Bus-13	0.0215	0.0707
Bus-52	Bus-54	0.0744	0.2444
Bus-52	Bus-56	0.0584	0.0349
Bus-53	Bus-56	0.0212	0.0834
Bus-77	Bus-83	0.0132	0.0437
Bus-83	Bus-88	0.0454	0.1801
Bus-85	Bus-87	0.0123	0.0505
Bus-85	Bus-90	0.0252	0.1170
Bus-9	Bus-30	0.0474	0.01563
Bus-26	Bus-38	0.0108	0.0331
Bus-38	Bus-68	0.0180	0.0919
Bus-64	Bus-68	0.0258	0.1170

Table A2

Detailed critical loads information.

At Bus	Total Load (MW)	Critical Load (MW)	Critical Level	χ L
Bus-1	51	20	4	1200
Bus-3	39	15	1	300
Bus-4	39	10	2	600
Bus-7	19	10	2	600
Bus-11	70	30	2	600
Bus-15	90	42	1	300
Bus-16	25	10	2	600
Bus-19	45	30	2	600
Bus-27	71	35	2	600
Bus-29	24	10	2	600
Bus-31	43	28	3	900
Bus-33	23	20	3	900
Bus-36	31	15	2	600
Bus-40	66	10	2	600
Bus-42	96	15	2	600
Bus-43	18	10	2	600
Bus-45	53	8	2	600
Bus-46	28	10	2	600
Bus-47	34	15	1	300
Bus-49	87	40	1	300
Bus-54	113	80	2	600
Bus-55	22	40	4	1200
Bus-56	18	50	1	300
Bus-59	277	20	2	600
Bus-60	78	40	2	600
Bus-62	77	40	2	600
Bus-66	39	10	3	900
Bus-70	66	40	2	600
Bus-74	68	50	3	900
Bus-75	47	20	2	600
Bus-76	68	50	2	600
Bus-77	61	60	4	1200
Bus-78	71	20	2	600
Bus-80	130	90	2	600
Bus-82	54	20	3	900
Bus-85	24	10	2	600
Bus-88	48	30	2	600
Bus-90	163	50	3	900
Bus-92	65	15	2	600
Bus-95	42	20	2	600
Bus-98	34	10	2	600
Bus-100	37	35	4	1200
Bus-101	22	10	2	600
Bus-103	23	10	2	600
Bus-104	38	15	1	300
Bus-105	31	10	1	300
Bus-107	50	14	1	300
Bus-110	39	20	1	300
Bus-115	22	10	2	600

Table A3

The specific information of cable lines.

Cable Lines
From Bus	To Bus
Bus-3	Bus-11
Bus-4	Bus-11
Bus-15	Bus-19
Bus-29	Bus-31
Bus-30	Bus-38
Bus-45	Bus-49
Bus-54	Bus-55
Bus-63	Bus-64
Bus-77	Bus-80
Bus-100	Bus-104
Bus-105	Bus-107
Bus-109	Bus-110

References

1. Alexandre, M.; Goran, S.; Rodrigo, M.; Alexandre, S.; Loannis, K. A Five-Level MILP Model for Flexible Transmissio-n Network Planning Under Uncertainty: A Min–Max Regret Approach. IEEE Trans. Power Syst.; 2017; 33, pp. 486-501. [DOI: https://dx.doi.org/10.1109/TPWRS.2017.2710637]

2. Alexandre, M.; David, P.; Alexandre, S.; Enzo, S. Reliable Renewable Generation and Transmission Expansion Plannin-g: Co-Optimizing System’s Resources for Meeting Renewable Targets. IEEE Trans. Power Syst.; 2017; 32, pp. 3246-3257. [DOI: https://dx.doi.org/10.1109/TPWRS.2016.2631450]

3. Jia, L.; Zuyi, L.; Feng, L.; Hongxing, Y.; Xuemin, Z.; Shengwei, M.; Naichao, C. Robust Coordinated Transmission an-d Generation Expansion Planning Considering Ramping Requirements and Construction Periods. IEEE Trans. Power Syst.; 2018; 36, pp. 268-280. [DOI: https://dx.doi.org/10.1109/TPWRS.2017.2687318]

4. Weixin, Z.; Changzheng, S.; Bo, H.; Kaigui, X.; Pierluigi, S.; Maosen, C. Transmission Defense Hardening Against Typ-hoon Disasters Under Decision-Dependent Uncertainty. IEEE Trans. Power Syst.; 2023; 38, pp. 1653-2665. [DOI: https://dx.doi.org/10.1109/TPWRS.2022.3194307]

5. Tao, D.; Ming, Q.; Zekai, W.; Bo, C.; Chen, C.; Mohammad, S. Power System Resilience Enhancement in Typhoons Using a Three-Stage Day-Ahead Unit Commitment. IEEE Trans. Power Syst.; 2021; 12, pp. 2153-2164. [DOI: https://dx.doi.org/10.1109/TSG.2020.3048234]

6. Mahdavi, M.; Antunez, C.S.; Ajalli, M.; Romero, R. Transmission Expansion Planning: Literature Review and Classific-ation. IEEE Syst. J.; 2019; 13, pp. 3129-3140. [DOI: https://dx.doi.org/10.1109/JSYST.2018.2871793]

7. Zhenzhi, L.; Fushuan, W.; Huifang, W.; Guanqiang, L.; Xiaojun, Y. CRITIC-Based Bus Importance Evaluation in Skele-ton-Network Reconfiguration of Power Grids. IEEE Trans. Circuits Syst. II Exp. Briefs.; 2018; 65, pp. 206-210. [DOI: https://dx.doi.org/10.1109/TCSII.2017.2703989]

8. Zhang, H.; Heydt, G.T.; Vittal, V.; Quintero, J. An improved network model for transmission expansion planning consi-dering reactive power and network losses. IEEE Trans. Power Syst.; 2013; 28, pp. 3471-3479. [DOI: https://dx.doi.org/10.1109/TPWRS.2013.2250318]

9. Yilin, X.; Ying, X. Transmission Expansion Planning Considering Wind Power and Load Uncertainties. Energies; 2022; 15, 7140. [DOI: https://dx.doi.org/10.3390/en15197140]

10. Yan, L.; Xueping, G. Skeleton-Network Reconfiguration Based on Topological Characteristics of Scale-Free Networks a-nd Discrete Particle Swarm Optimization. IEEE Trans. Power Syst.; 2007; 22, pp. 1267-1274. [DOI: https://dx.doi.org/10.1109/TPWRS.2007.901486]

11. Gu, X.; Zhong, H. Optimization of Network Reconfiguration Based on A Two-Layer Unit-Restarting Framework for P-ower System Restoration. IET Gener. Transm. Distrib.; 2012; 6, pp. 693-700. [DOI: https://dx.doi.org/10.1049/iet-gtd.2011.0591]

12. Zhenzhi, L.; Fushuan, W.; Yusheng, X. A Restorative Self-Healing Algorithm for Transmission Systems Based on Complex Network Theory. IEEE Trans. Smart Grid.; 2016; 7, pp. 2154-2162. [DOI: https://dx.doi.org/10.1109/TSG.2016.2539199]

13. Joshua, A.T.; Franz, S.H. Linear Relaxations for Transmission System Planning. IEEE Trans. Power Syst.; 2011; 26, pp. 2533-2538. [DOI: https://dx.doi.org/10.1109/TPWRS.2011.2145395]

14. Romero, R.; Asada, E.N.; Carreno, E.; Rocha, C. Constructive Heuristic Algorithm in Branch-And-Bound Structure Ap-plied to Transmission Network Expansion Planning. IET Gen. Trans. Distrib.; 2007; 1, pp. 318-323. [DOI: https://dx.doi.org/10.1049/iet-gtd:20060239]

15. Verma, A.; Panigrahi, B.K.; Bijwe, P.R. Harmony Search Algorithm for Transmission Network Expansion Planning. IET Gen. Trans. Distrib.; 2010; 1, pp. 663-673. [DOI: https://dx.doi.org/10.1049/iet-gtd.2009.0611]

16. Leite da Silva, A.M.; Rezende, L.S.; Honorio, L.M.; Manso, L.A.F. Performance Comparison of Metaheuristics to Solv-e The Multi-Stage Transmission Expansion Planning Problem. IET Gen. Trans. Distrib.; 2011; 1, pp. 360-367. [DOI: https://dx.doi.org/10.1049/iet-gtd.2010.0497]

17. Neeraj, G.; Mahdi, K.; Nilesh, P.; Tomanobu, S. A Bi-Level Evolutionary Optimization for Coordinated Transmission Expansion Planning. IEEE Access.; 2018; 28, pp. 48455-48477. [DOI: https://dx.doi.org/10.1109/ACCESS.2018.2867954]

18. Mohsen, R.; Ruben, R.; Marcos, J.R. Strategies to Reduce the Number of Variables and the Combinatorial Search Space of the Multistage Transmission Expansion Planning Problem. IEEE Trans. Power Syst.; 2013; 28, pp. 2164-2173. [DOI: https://dx.doi.org/10.1109/TPWRS.2012.2223241]

19. Jin, Y.X.; Cheng, H.Z.; Yan, H.Y.; Zhang, L. New Discrete Method for Particle Swarm Optimization and Its Applicatio-n in Transmission Network Expansion Planning. Elect. Power Syst. Res.; 2007; 77, pp. 227-233. [DOI: https://dx.doi.org/10.1016/j.epsr.2006.02.016]

20. Pouria, M.; Seyed, H.H.; Majid, O.B.; Mohammad, S. A Scenario-Based Multi-Objective Model for Multi-Stage Transmis-sion Expansion Planning. IEEE Trans. Power Syst.; 2011; 26, pp. 470-478. [DOI: https://dx.doi.org/10.1109/TPWRS.2010.2048930]

21. Dong, W.; Chen, X.; Yang, Q. Data-driven scenario generation of renewable energy production based on controllable g-enerative adversarial networks with interpretability. Appl. Energy; 2022; 308, 118387. [DOI: https://dx.doi.org/10.1016/j.apenergy.2021.118387]

22. Zhan, J.; Chung, C.; Zare, A. A fast solution method for stochastic transmission expansion planning. IEEE Trans. Power Syst.; 2017; 32, pp. 4684-4695. [DOI: https://dx.doi.org/10.1109/TPWRS.2017.2665695]

23. Ma, Z.; Gao, J.; Hu, W.; Dinavahi, V. Risk-adjustable stochastic schedule based on sobol augmented latin hypercube sa-mpling considering correlation of wind power uncertainties. IET Renew. Power Gener.; 2021; 15, pp. 2356-2367. [DOI: https://dx.doi.org/10.1049/rpg2.12169]

24. Bukenberger, J.P.; Webster, M.D. Approximate latent factor algorithm for scenario selection and weighting in transmissi-on expansion planning. IEEE Trans. Power Syst.; 2020; 35, pp. 1099-1108. [DOI: https://dx.doi.org/10.1109/TPWRS.2019.2942925]

25. Yu, H.; Chung, C.Y.; Wong, K.P.; Zhang, J.H. A Chance Constrained Transmission Network Expansion Planning Metho-d with Consideration of Load and Wind Farm Uncertainties. IEEE Trans. Power Syst.; 2009; 24, pp. 1568-1576. [DOI: https://dx.doi.org/10.1109/TPWRS.2009.2021202]

26. Deping, K.; Chung, C.Y.; Yuanzhang, S. A Novel Probabilistic Optimal Power Flow Model with Uncertain Wind Pow-er Generation Described by Customized Gaussian Mixture Model. IEEE Trans. Sustain. Energy.; 2016; 7, pp. 200-212. [DOI: https://dx.doi.org/10.1109/TSTE.2015.2489201]

27. Valverde, G.; Saric, A.T.; Terzija, V. Probabilistic Load Flow with Non-Gaussian Correlated Random Variables Using G-aussian Mixture Models. IET Gen. Trans. Distrib.; 2012; 6, pp. 701-709. [DOI: https://dx.doi.org/10.1049/iet-gtd.2011.0545]

28. Chen, J.F.; Sun, X.; Duan, X.Z.; Yang, Z.L.; Zhou, H.B.; Wang, J.B. A Chance-constrained Approach for Available Tra-nsfer Capability Evaluation for Power Systems with Wind Farm Integration. Proc. CSEE.; 2019; 39, pp. 6804-6814. [DOI: https://dx.doi.org/10.13334/j.0258-8013.pcsee.182088]

29. Hong, H.F.; Hu, Z.S.; Guo, R.P.; Ma, J.; Tian, J. Directed Graph-Based Distribution Network Reconfiguration for Oper-ation Mode Adjustment and Service Restoration Considering Distributed Generation. J. Mod. Power Syst. Clean. Energy.; 2017; 5, pp. 142-149. [DOI: https://dx.doi.org/10.1007/s40565-016-0198-3]

30. Baldwin, T.L.; Mili, L.; Phadke, A.G. Dynamic Ward Equivalents for Transient Stability Analysis. IEEE Trans. Power Syst.; 1994; 9, pp. 59-67. [DOI: https://dx.doi.org/10.1109/59.317557]

31. Shuwei, X.; Wenchuan, W.; Tao, Z.; Zhenyi, W. Convex Relaxation Based Iterative Solution Method for Stochastic Dynamic Economic Dispatch with Chance Constrain. Autom. Elect. Power Syst.; 2020; 44, pp. 43-51. [DOI: https://dx.doi.org/10.7500/AEPS20190906001]

32. Ahmed, S. Convex Relaxations of Chance Constrained Optimization Problems. Optim. Letter.; 2014; 8, pp. 1-2. [DOI: https://dx.doi.org/10.1007/s11590-013-0624-7]

33. Elia Group: Wind Power Generation. Available online: https://www.elia.be/en/grid-data/power-generation/wind-power-generation (accessed on 28 June 2023).