Renewable Generation and Transmission Expansion Planning Coordination with Energy Storage System: A Flexibility Point of View

: This paper presents a method for coordinated network expansion planning (CNEP) in which the difference between the total cost and the ﬂexibility beneﬁt is minimized. In the proposed method, the generation expansion planning (GEP) of wind farms is coordinated with the transmission expansion planning (TEP) problem by using energy storage systems (ESSs) to improve network ﬂexibility. To consider the impact of the reactive power in the CNEP problem, the AC power ﬂow model is used. The CNEP constraints include the AC power ﬂow equations, planning constraints of the different equipment, and the system operating limits. Therefore, this model imposes hard nonlinearity onto the problem, which is linearized by the use of ﬁrst-order Taylor’s series and the big-M method as well as the linearization of the circular plane. The uncertainty of loads, the energy price, and the wind farm generation are modeled by scenario-based stochastic programming (SBSP). To determine the effectiveness of the proposed solution approach, it is tested on the IEEE 6-bus and 24-bus test systems using GAMS software.


Motivation and Approach
In recent years, the expansion planning problem of power systems has been faced with new challenges. In the conventional expansion planning models, the expansion planning of transmission systems and the generation expansion planning were two separate problems [1,2]. Renewable energy sources (RESs) are becoming increasingly popular due to their lower operating costs and lower environmental impact. Thus, determining the optimal size and location of RESs in the planning of power systems is important [3,4]. However, there are significant challenges associated with the operation of RESs. Due to their uncertain nature, they have negative impacts on the flexibility of the system [5]. The flexibility term is defined as "the modification of generation or/and consumption patterns in reaction to an external price or an activation signal to provide a service in the electrical system" [5]. One of the main solutions to overcome the inflexibility problem of a power system is the deployment of flexible sources (FSs) such as energy storage systems (ESSs) [6]. Therefore, the expansion planning of transmission systems can no longer be considered a separate problem from generation expansion planning, due to the system's flexibility challenges.

Contributions
Based on the previous research works in this area, the research gaps can be summarized as follows:

•
None of the previous studies on CNEP consider the flexibility index of power systems. However, it is important in a power system in the presence of RESs to modify the generation injection and/or consumption patterns that are not suitable from a power balancing viewpoint due to the uncertainties of RESs. This index can be improved in the power system using the optimal planning of energy storage systems such as batteries. • Many research works reviewed above take into account DC power flow constraints in their CNEP models. This simplification results in the elimination of reactive power studies in an expansion planning problem. It should be noted that the capacity of lines and power sources also depends on the reactive power as well as the active power, and without considering the reactive power, the CNEP solution may lead to impractical results. • Some works use evolutionary methods to solve the CNEP problem. These methods search for different directions of the solution space, and their progress is usually slow. Therefore, the computation time of evolutionary methods is high. In addition, their convergence to the global optimum for the CNEP problem cannot be guaranteed.
To address the above issues, this paper presents a CNEP method to coordinate the optimal expansion planning of transmission systems with the size and the placement of flexible and renewable sources, i.e., ESSs and wind farms, considering the network's flexibility. In the proposed model, the difference between investment and operation costs and the flexibility benefit is minimized. The linearized AC optimal power flow equations, ESS and RES operational constraints, and the flexibility limit are the constraints of the proposed model. Hence, this paper presents a mixed-integer linear programming (MILP) model to guarantee the global optimal solution with a low computational burden. In addition, the uncertain parameters of the problem are modeled based on the scenariobased stochastic programming (SBSP) approach. Therefore, the roulette wheel mechanism (RWM) is implemented for scenario generation, and the Kantorovich method is used for the scenario reduction to model the forecasted error uncertainty of the load, the energy price, and the RES power. The main contributions of this paper can be summarized as follows:

•
Modeling the coordinated network expansion planning in the power system according to a hybrid method of optimal planning of renewable and flexible sources and TEP to improve the operation and the flexibility indices. • Obtaining the linear model of the proposed strategy, considering the least calculation error concerning the original method using linearized AC optimal power flow equations. • Using scenario-based stochastic programming based on a hybrid approach with the roulette wheel mechanism and Kantorovich method to model the uncertainty of the load, the energy price, and the RES power.

Paper Organization
The rest of this paper is organized as follows: the model of stochastic MINLP and the MILP of the CNEP strategy are presented in Section 2. In Section 3, numerical results of the application of the proposed problem on test networks are shown, and finally, some concluding remarks are made in Section 4.

Stochastic CNEP Problem
In this section, the original mathematical model of the CNEP problem based on Figure 1 is presented. In the proposed CNEP model, the investment and operation costs of power systems and the benefit gained from increasing the flexibility are included. In the following subsections, the mathematical formulations are discussed.

Stochastic CNEP Problem
In this section, the original mathematical model of the CNEP problem based on Figure 1 is presented. In the proposed CNEP model, the investment and operation costs of power systems and the benefit gained from increasing the flexibility are included. In the following subsections, the mathematical formulations are discussed.  (A) Objective function: the objective function of the CNEP problem is expressed in Equation (1), and consists of five terms. The first one indicates the operation cost of generation units, and the second one presents the annual investment cost of the transmission expansion. The third and fourth terms refer to the RESs' annual investment cost, and annual investment and operation cost of ESSs, respectively, and the flexibility index is formulated in the last term of the objective function. An annual step time in the horizon time is considered in this paper. Thus, the operation cost is multiplied by 8760, i.e., 24 × 365. It should be noted that the first and fifth terms of the objective function and the operation cost of ESS depend on uncertain parameters, while the second term, third term, and ESS investment costs do not. This means that the generating and ESS schedules are associated with the scenarios' realization, as well as the system flexibility; however, the investment in new transmission lines, ESS, and RES cannot be changed by each scenario [32]. (A) Objective function: the objective function of the CNEP problem is expressed in Equation (1), and consists of five terms. The first one indicates the operation cost of generation units, and the second one presents the annual investment cost of the transmission expansion. The third and fourth terms refer to the RESs' annual investment cost, and annual investment and operation cost of ESSs, respectively, and the flexibility index is formulated in the last term of the objective function. An annual step time in the horizon time is considered in this paper. Thus, the operation cost is multiplied by 8760, i.e., 24 × 365. It should be noted that the first and fifth terms of the objective function and the operation cost of ESS depend on uncertain parameters, while the second term, third term, and ESS investment costs do not. This means that the generating and ESS schedules are associated with the scenarios' realization, as well as the system flexibility; however, the investment in new transmission lines, ESS, and RES cannot be changed by each scenario [32].
(B) AC power flow equations: These equations are introduced in (2)-(5), which correspond to the active and reactive power balance in each bus and time, (2) and (3), and the active and reactive power flow of lines in each time, (4) and (5), respectively. It should be noted that the multiplication of binary variable x in (4) and (5) is due to modeling the effect of the line construction. These constraints should be exerted only for the line candidates that are selected to be constructed. Also, note that the line between buses i and j is the same as the one between buses j and i. This fact is modeled by (6). In addition, A is a binary matrix whose elements represent the existence of transmission lines between different busses [21].
PL i,j,t,y,ω = g i,j V i,t,y,ω 2 − V i,t,y,ω V j,t,y,ω g i,j cos(θ i,t,y,ω − θ j,t,y,ω ) (C) System operation limitations: The operational limits of the system, i.e., the line capacity, the generating unit capacity, and the voltage limits, are expressed through (7)-(9) by the power flow equations [6].
(D) RES constraints: Constraint (10) presents the RES active power, where the maximum available power of RES (i.e., PR max ) is the function of wind speed/solar radiation for a wind/solar farm. Additionally, in constraint (10), if the RES i is constructed, (E) Flexibility resource constraints: The flexibility sources are ESSs in this paper due to their fast response characteristics [28], and hence, constraints (11)-(16) present the ESS operational equations. These constraints are charge rate limit, (11); discharge rate limit, (12); logical limit avoiding simultaneous ESS charging and discharging operation, (13); energy stored in an ESS, (14); minimum and maximum energy limits, (15); and the initial energy of an ESS, (16). Additionally, in constraints (13), (15), and (16) (F) System flexibility constraints: This paper uses a metric to quantify the technical flexibility level of both individual flexible resources and the whole system based on the approach in [29]. This metric quantifies the technical ability of a system or resources to cope Appl. Sci. 2021, 11, 3303 6 of 18 with the flexibility requirement resulting from the variability and the uncertainty produced by RESs. Flexible resources, i.e., ESSs, can be used to meet the system flexibility requirement. Also, there are upward and downward flexibility indices for each resource [29]. If the output power of flexible resources in the scenario ω is greater than its output power in the base case, there is upward flexibility for it. However, there is downward flexibility for the flexible resource if its output power in the scenario ω is less than the output power of this source in the base case [29]. The upward and downward flexibility for ESS can be obtained from (17). Moreover, system flexibility (SF) is calculated based on (18) [29].

The Linear CNEP Model
As can be observed from the aforementioned problem model, constraints (4), (5), (7) and (8) are nonlinear. In addition, Equations (4) and (5) are non-convex due to the terms of cos(θ i,t,y,ω − θ j,t,y,ω ) and sin(θ i,t,y,ω − θ j,t,y,ω ) [33][34][35], and they include a binary variable. Therefore, the model of (1)-(18) is a non-convex mixed-integer nonlinear problem (MINLP). This will slow down the solving procedure and may result in infeasible solutions in actual large-scale networks [36][37][38]. In this section, an equivalent linear model is proposed to speed up the solving procedure and guarantee the global optimality of the solution.
(A) The Linear approximation of AC power flow constraints: Nonlinear constraints in this section are (4) and (5). For linearizing these equations, the following assumptions are considered [38,39]:

•
The difference in the voltage phases of the two ends of a transmission line is under 6 degrees or 0.105 radians.

•
Since the magnitude of bus voltages in the transmission network should be between 0.95 and 1.05 per unit, the magnitude of bus voltages is very close to 1 per unit.
According to the mentioned assumptions, the bus voltages can be expressed as 1 + ∆V where the value of ∆V is much smaller than 1. Also, based on the first assumption, cos(θ i,t,y,ω − θ j,t,y,ω ) and sin(θ i,t,y,ω − θ j,t,y,ω ) are equal to 1 and (θ i,t,y,ω − θ j,t,y,ω ), respectively. It is worth noting that the values of ∆V 2 , ∆V × (θ i,t,y,ω − θ j,t,y,ω ) and ∆V i × ∆V j are small enough to ignore, which is the case in this paper. Hence, Equations (4) and (5) are rewritten as below.
In order to linearize the multiplication of x i,j with other continuous terms, the big-M method can be used here as follows: ∀i, j, t, y, ω where M is a large enough number which is assumed to be 10 6 . Also, the linearization procedure of Equation (23) is described in the next section.
(B) The linear approximation of system operation limitations: Equations (7), (8) and (23) are nonlinear. These relations represent circular inequalities, which denote the variation interval of active and reactive power inside a circle with a radius equal to the maximum apparent power. Nonetheless, determination of the linear equations corresponding to the circular plane with Equation (24) can be performed as illustrated in Figure 2 [36][37][38]. Hence, the circular plane is obtained by unifying several square planes. For example, based on Figure 2, the linear description of the circular plane is as in Equation (25).
(B) The linear approximation of system operation limitations: Equations (7), (8) and (23) are nonlinear. These relations represent circular inequalities, which denote the variation interval of active and reactive power inside a circle with a radius equal to the maximum apparent power. Nonetheless, determination of the linear equations corresponding to the circular plane with Equation (24) can be performed as illustrated in Figure 2 [36][37][38]. Hence, the circular plane is obtained by unifying several square planes. For example, based on Figure 2, the linear description of the circular plane is as in Equation (25).  According to Figure 2, Equation (24) has a large linearization error. Thereby, to reduce the linearization error, the number of square planes with different angles from the horizontal axis should be increased. Thus, at first 360 degrees of circle perimeter are divided into equal parts ∆α. Then, the line equation is linearized for each k × ∆α, where k is the counter of the linearized part. Finally, the calculated line equation is integrated with a circle with a radius less than or equal to r. Therefore, the linear approximation of the circular plane is expressed as below: x cos(k∆α) + y sin(k∆α) ≤ r k ∈ ϕ k = {0, 1, . . . , n k − 1} (26) where n k denotes the number of linearized parts. For example, if 180 square planes are used for the linearization of a circular plane, n k and ∆α will be equal to 180 and 2, respectively. Consequently, the planes x ≤ r and y ≤ r correspond to k = 0 and k = 45. Therefore, the linear approximation of (23) and (8) is as follows.
Subject to: As a final remark, in order to transform the model in (29)-(31) into a CNEP problem with DC power flow equations, the terms Q, ∆V, and g should be set to zero, and Equations (3), (22) and (31) must be removed. The term cos(k × ∆α) should be deleted from equations.

Uncertainties Model
In the proposed problem, the load parameters, i.e., PD and QD, the charging price (λ), and the maximum RES active power (PR max ) are considered to be the sources of uncertainty. In this paper, the roulette wheel mechanism (RWM) [27] for the generation of scenario samples is used for this case. It is noted that the RWM generates scenarios for load and price forecasting errors based on the normal distribution function [39], for the maximum RES active power based on Weibull distribution functions [27]. The Kantorovich method [27] is selected in this paper to reduce the number of scenarios, where the proposed scenario generation/reduction algorithm details are presented in [27].

Numerical Results and Discussion
To show the effectiveness of the proposed CNEP problem with AC power flow equations, it is tested on the IEEE 6-bus and IEEE 24-bus test systems, which are usually considered as test cases in previous research works in this area, such as [40][41][42]. Additionally, the proposed method is implemented in the GAMS software package using the solver of CPLEX/BONMIN for the MILP/MINLP problem [43]. An explanation about these solvers can be found in [44][45][46]. Note that the total number of linearization segments of circular constraint is considered to be 180 to obtain a low calculation error.

IEEE 6-Bus Test System
The diagram of the IEEE 6-bus test system is shown in Figure 3 [47], where bold lines denote the existing lines, and dashed lines represent the proposed lines to be constructed in the network. The power demand at the peak load hour, line data, and generator characteristics are described in Tables 1-3. Also, a planning horizon of 5 years is considered for the transmission lines. Table 4 represents the load percentage for each year, and the coincidence factor is considered to be 0.7 [21]. Moreover, it is forecasted that the wind farm with a capacity of 50 MW, and ESS with 200 MWh and charging/discharging rate of 50 MW and 88% efficiency can be installed in bus 4. Additionally, the forecasted daily load factor, the energy price, and the maximum RES (wind farm) power percentage are plotted in Figure 4 [34,48]. rate of 50 MW and 88% efficiency can be installed in bus 4. Additionally, the forecasted daily load factor, the energy price, and the maximum RES (wind farm) power percentage are plotted in Figure 4 [34,48].
The base power and base voltage are equal to 100 MW and 230 kV, and the minimum and maximum voltage is 0.95 and 1.05 pu. Moreover, the reference bus is bus 1 with a voltage of 1.05∠0 pu.   Table 2. Line characteristics of 6-bus network [47].        In this section, the deterministic TEP model considers only peak load hours, and both ESS and RES are removed from the proposed problem. The corresponding results for the three models are presented in Figure 5 and Table 5. In this section, the deterministic TEP model considers only peak load hours, and both ESS and RES are removed from the proposed problem. The corresponding results for the three models are presented in Figure 5 and Table 5.

Capacity (MVA) Construction Cost (M$) Reactance (pu) Resistance (pu) From-To Bus
According to this figure, the results of modeling of the TEP problem based on the DC power flow, the nonlinear power flow, and the linear approximation of AC power flow are the same in selecting the transmission lines for construction. For each of the three models, a transmission line is constructed between busses 4 and 6 by considering the load increase. However, it is seen that in Figure 5, the active and reactive power of generators and bus voltages have different values for different TEP models. The difference in bus voltages is negligible where linear or nonlinear AC power flow is used, but the difference in generators' active and reactive power is about 3%. Also, it is observed that in the model based on DC power flow, bus voltages are always equal to 1 per unit and the reactive power is always zero.     Table 5 represents the obtained results for the operation and investment cost, constructed lines, network active and reactive losses, computation time, and the situation of problem-solving for three models. Based on this table, the construction cost is equal for the three models, due to their unique configuration. The exploitation cost and network's active and reactive losses during the planning are equal in DC and linear AC power flow (models 1 and 3). Based on the obtained results, network losses in these two models are zero, due to the equations that are used in them. Since modeling based on nonlinear AC power flow (model 2) considers network losses, the operation cost of this model is greater than the two other ones. However, it is worth noting that model 2 as an MINLP problem achieves more accurate results, but its computation time is significantly larger than the computation time of models 1 and 3. Based on Table 6, model 2 does not have a globally optimal solution. Noticeably, the proposed TEP problem using linear AC power flow equations (model 3) possesses a slight approximation in comparison with model 2. It has a very small computation time and achieves a globally optimal solution. Also, the computation time of model 1 is short and it obtains an optimal solution, but the obtained result of this model is not feasible because the reactive power is not considered in DC power flow equations. Computation time (s) 5

Solution situation Optimal
Herein, the results of the stochastic TEP method with the standard deviation of 10% for active and reactive load uncertainty are shown in Figure 6 and Table 6. By comparing Figures 5b and 6, it is observed that in the stochastic model, two transmission lines are proposed to be constructed which include 2-6 and 4-6. Also, the expected active and reactive power of generators has been increased in comparison with the deterministic model, and the expected load bus voltages have been reduced. Based on Table 6, the total cost of operation and construction is about 153.368 M$, which will be equal to 93.062 M$ in the deterministic problem. The cost difference is 60.306 M$, which can be considered the prediction cost or the uncertainty cost or the security cost due to considering the load uncertainty.    In this paper, the term FC is considered 2000 $ to obtain high flexibility in a 6-bus network. The results of three cases are presented in Table 7, and accordingly, it is seen that the proposed CNEP strategy has the following results:

•
Reduce the operation cost: The operation cost in case II is 3,045,670 $, a reduction of 9.57% with respect to the TEP model. This term is 13.13% for the coordinated renewable sources and transmission expansion planning method.

•
Reduce the investment cost: In the proposed CNEP method (case III), there is one constructed line that is between busses 2 and 6, while the TEP model includes two constructed lines. • Improve the voltage profile: In case III, the maximum voltage drop is 0.053 p.u., while it is 0.065 and 0.059 p.u. in cases I and II, respectively. Therefore, the CNEP strategy can improve the voltage profile by about 18.46%.

•
Obtain high flexibility: Based on Table 7, the CNEP method includes system flexibility of 10.93 as well as 0.2186 million $ flexibility benefit, where it is high concerning cases I and II. In this paper, the term FC is considered 2000 $ to obtain high flexibility in a 6-bus network. The results of three cases are presented in Table 7, and accordingly, it is seen that the proposed CNEP strategy has the following results:

•
Reduce the operation cost: The operation cost in case II is 3,045,670 $, a reduction of 9.57% with respect to the TEP model. This term is 13.13% for the coordinated renewable sources and transmission expansion planning method.

•
Reduce the investment cost: In the proposed CNEP method (case III), there is one constructed line that is between busses 2 and 6, while the TEP model includes two constructed lines. • Improve the voltage profile: In case III, the maximum voltage drop is 0.053 p.u., while it is 0.065 and 0.059 p.u. in cases I and II, respectively. Therefore, the CNEP strategy can improve the voltage profile by about 18.46%. • Obtain high flexibility: Based on Table 7, the CNEP method includes system flexibility of 10.93 as well as 0.2186 million $ flexibility benefit, where it is high concerning cases I and II.

IEEE 24-Bus Network
The 24-bus test system structure is presented in [14]. This network has 34 transmission lines, 24 of which exist, and 10 others of which are proposed for construction. The proposed lines for construction are shown in Table 8. Moreover, the wind farms with a capacity of 150 MW and ESSs with the maximum capacity of 400 MWh, 100 MW charging/discharging rate, and 88% efficiency can be connected to busses 13, 17, 21, and 23. The other data such as load factor, energy price, maximum RES power, etc., are the same as 6-bus test network data. The base power and base voltage of this network are 100 MVA and 230 kV, respectively. The results for the proposed CNEP method are shown in Figures 7 and 8. Figure 7 shows the energy profile of ESSs, in which the charging/discharging operation mode is represented by increasing/decreasing ESS energy. Accordingly, it is charged at the period of 1:00 to 8:00 due to the low energy price, and it is discharged at the period of 17:00 to 22:00 due to the high energy price as shown in Figure 4b. In other words, ESSs are charged during the low energy price period and are discharged at the high energy price time to have a low operation cost, based on objective function (1). Moreover, Figure 8 expresses the expected flexibility benefit and expected operation cost for the IEEE 24-bus network. According to Figure 8a, the expected flexibility benefit is increased by increasing the flexibility price (FC) to 10 5 $, and it is constant to 2.5 × 10 5 if FC changes between 10 5 and 10 6 $. Also, the expected operation cost will be increased if FC increases between 0 and 10 5 $, and it is constant to 60 M$ if FC changes between 10 5 and 10 6 $. It is noted that the increase in FC causes the number of flexible sources or ESSs to be increased in the power system to obtain the high flexibility benefit; hence, the operation cost will be increased because a high number of ESSs is used in the power system in this condition. As a result, the CNEP can provide a transmission network with a lot of flexibility by setting a high flexibility incentive price. A suitable cost of operation value is also required. Finally, Table 9 presents the values of maximum voltage drop, system flexibility, and problem convergence. Accordingly, the maximum voltage drop from 1.05 p.u. is 0.064 p.u., which is less than the allowed value, i.e., 0.1 (1.05-0.95) p.u. Also, system flexibility is 19.4. In the end, these capabilities are calculated in 21 s, where the solution situation is optimal. Thus, the proposed method can obtain a reliable optimal solution in a low calculation time.

Conclusions
This paper presents a method for the CNEP problem that minimizes the difference between investment and operation costs and flexibility benefit according to optimal power flow equations and planning constraints of different equipment. At first, the linear formulation of the AC optimal power flow model is used in the proposed strategy. Stochastic programming according to RWM-based scenario generation and Kantorovich method-based scenario reduction is used in this paper to model the forecasted error uncertainty of the load, energy price, and RES. Numerical results show that the proposed configuration structure by using MINLP and MILP models with linear AC and DC power flows is the same. The difference between the three models is in the network parameters, such as bus voltages and losses, the produced active and reactive power, and the operation cost. This difference is caused by the zero loss in MILP with AC power flow, and also the zero loss and reactive power and the 1 p.u. voltage in the MILP model with DC power flow. The difference in the first model is lower than in the second one. Also, two of the models have very short computation time and optimal solutions, while the MINLP model has a very long computation time and its solution is not an optimal one. The mentioned items indicate the superior performance of the MILP model by using linear AC power flow compared with the two other models. Thereby, the proposed CNEP strategy can improve the network indices such as the voltage profile, and it can allow high flexibility in the test system.

Conflicts of Interest:
The authors declare no conflict of interest.

ch, dch
Binary variable of charging and discharging state inflexible source (without unit) E Stored energy in the flexible source F + , F − Upward and downward flexibility active power PF ch , PF dch Charging and discharging active power of flexible source PG, QG Active and reactive power of generation unit PL, QL Active and reactive power flow of line PR Active power of RES SF System flexibility (without unit) V, ∆V, θ Magnitude, deviation, and angle (in radian) of voltage x, x R , x F The binary variable of line, RES, and flexible source construction (without unit) (2) Constants: A Bus incidence matrix (if a line existed between buses i and j, A b,j is equal to 1, or otherwise zero) C Energy or operation price in $/MWh CF Coincidence factor E ini Initial energy of flexible source in p.u. E min , E max Minimum and maximum energy of flexible source in p.u.

FC
Flexibility price ($) g, b Conductance and susceptance of a line in p.u. IC L , IC R , IC F Investment cost of line, RES, and flexible source in $ PD, QD Active and reactive load in p.u.

PF max
Maximum charge/discharge rate in p.u.

PR max
Maximum RES active power in p.u. V min , V max Minimum and maximum value of voltage magnitude in p.u.

SL max
Maximum capacity of line in p.u.

SG max
Maximum capacity of generation unit in p.u.