In recent years, the incorporation of renewable energy (RE) has developed rapidly in power systems due to RE’s environmentally friendly characteristics. RE can help reduce carbon emissions from power systems; thus, RE penetration has become an important target in power systems [1]. For example, the European Union has issued the European Green Deal to increase RE penetration. However, the RE from wind and solar power is affected by climate change and has strong randomness and volatility. When large-scale RE is integrated into the network, its uncertainties will have a significant impact on the balance of the power system. As a long-term decision-making method, power system planning is the primary strategy for coping with climate change in this context. It is only when renewable-dominant power systems obtain sufficient adjustable resources through planning that they will be able to effectively ensure a power balance.

As the key component of a power system, the transmission network can transmit RE power to the load side to ensure a power balance; thus, ensuring sufficient transmission capacity through transmission expansion planning (TEP) is necessary [[2], [3], [4]]. While TEP can ensure the integration and transmission of RE power, it cannot fully cope with the randomness and volatility of RE power. Therefore, components with flexible adjustment capabilities are still required to ensure a power balance. As a kind of common flexible resource, energy storage systems (ESSs) can alleviate transmission congestion and mitigate fluctuations in RE power, so they are expected to play an important role in power systems with high penetration of RE [5,6]. As a result, TEP problems involving EESs have received widespread attention. Ref. [7] investigates the problem of battery-based energy storage transportation in TEP to improve the utilization of RE resources, while Ref. [8] establishes a robust TEP model with ESSs to cope with the uncertainties of RE. Ref. [9] proposes a bi-level model to analyze the benefits of ESSs and transmission construction for RE stations, and Ref. [10] establishes a tri-level model to analyze investment strategies for ESSs and transmission lines in a market environment.

Based on the above analysis, coordinated transmission renewable–storage expansion planning (TRSEP) is an effective decision-making method to cope with climate change and achieve the RE development target. Power systems are generally not isolated and often interact with other power systems. When there is an electricity surplus or shortage in a power system, electricity trading between different power systems is conducive to ensuring a power system balance. For example, a power system with abundant power-generation resources can export excess electricity to other power systems through tie lines on the premise of ensuring the internal load demands of the system. Therefore, tie lines must be considered during planning as external loads of the export system.

In general, annual transaction contracts are signed between the export system and the import system. In the context of energy transformation, many power systems are placing greater emphasis on green power trading. The European Union has already adopted Guarantees of Origins to demonstrate the proportion of electricity supplied by RE to load demands. Moreover, some systems even require that the RE energy transmitted by tie lines meet a certain proportion. However, the literature has not considered the analysis of new energy transmission pathways in power system planning, leaving a research gap that is addressed in this paper.

In renewable-dominant power systems, analysis of the RE transmission path will receive increased attention, as it can provide a calculation basis for the settlement of green electricity trading. However, existing studies of power system planning with RE penetration targets typically consider the system’s internal loads without the external loads in tie lines. When considering external load demands, the RE penetration targets can differ for internal and external loads. In addition, there may be some areas within power systems that aim to receive more RE generation to meet their load demands. These areas, which can be referred to as “energy transition areas,” have higher requirements for the use of clean electricity than other areas within power systems. It is also necessary to clarify the RE transmission path in order to supply internal loads in power systems. Therefore, there is a need to calculate the RE transmission path in renewable-dominant power systems.

A power flow tracing method for the analysis of energy transmission paths is proposed for the first time in Ref. [11], where it is applied to calculate the shares of power sources in transmission facilities. In this method, the power flow follows the principle of proportional sharing. On this basis, Ref. [12] establishes a network-based carbon emission flow (CEF) model to analyze the CEF from generation to demand, which can be embedded in optimization problems. The CEF intensity (or density) is defined in order to calculate the carbon emissions of each branch and bus.

Research on CEF is relatively mature. Ref. [13] establishes a coordinated electricity and gas network planning model with CEF constraints, and Ref. [14] proposes low-carbon constrained multiple energy systems planning with CEF. Ref. [15] adopts the CEF model for the operation of multiple energy systems, while Ref. [16] applies CEF to formulate a scheduling plan for power plants in a market environment. Moreover, an intertemporal CEF model for ESSs is established in Ref. [17] that can describe the temporal variation of CEF in ESSs.

By referring to the definition of CEF, it is possible to similarly track the path of RE power. The RE power-flow density (RE-PFD) can be defined to reflect the proportion of RE power in the power of the electrical components. On this basis, the RE power transmitted to the internal and external loads in the TRSEP model can be calculated, thereby achieving both the RE development goals for the system and the requirements for green power trading between systems. Thus, this paper proposes a TRSEP model for renewable-dominant power systems that considers the transmission path of RE, as an effective decision-making method to cope with climate change.

Due to the introduction of the RE-PFD, the model may exhibit nonlinear terms, which will pose significant challenges for solving the model. In order to eliminate nonlinear terms, Ref. [18] adopts a piecewise linear approximation method, which discretizes the continuous feasible region by introducing a large number of 0–1 variables. Although this method eliminates nonlinear terms, it makes the model too complex to solve. Ref. [19] proposes a McCormick method that does not introduce any 0–1 variables, but it expands the feasible region and the obtained solution may not be feasible. This method is used in Ref. [20] to linearize the nonlinear terms that occur in power system planning, yielding an approximate solution. Based on Ref. [19], a segmented McCormick method was designed in Ref. [21], which reduces the feasible region to some extent; however, it also greatly increases the number of 0–1 variables and still cannot guarantee obtaining feasible solutions. Based on the above analysis, the McCormick method does not significantly increase the computational complexity, which is extremely advantageous for solving large-scale power system planning problems. If a correction strategy can be designed to enable the McCormick method to effectively obtain feasible solutions, the planning model proposed in this paper can be well solved.

In summary, current research does not consider the RE transmission path in power system planning models. In addition, in research on solving nonlinear optimization problems, existing methods may add 0–1 variables or fail to obtain feasible solutions. Therefore, research on methods that can effectively obtain feasible solutions is needed.

The specific contributions of this paper are as follows:

(1) By referring to the definition of CEF and extending the research in Ref. [11], the concept of RE-PFD is defined, which can provide a theoretical basis for the analysis of the RE transmission path.

(2) Based on the definition of RE-PFD, the relevant operational constraints of RE-PFD are proposed, which can describe the intertemporal characteristics of RE-PFD and clearly distinguish the RE transmission path into internal loads, external loads, and energy losses. Moreover, a complete TRSEP model to cope with the impacts of climate change and achieve the development target of RE penetration is established.

(3) To handle the bilinear terms generated by the power and RE-PFD, a customized feasibility correction strategy based on the McCormick method is designed. More specifically, to address the issue of the McCormick method expanding the feasible region and being unable to yield feasible solutions, a strategy is designed to reduce the feasible region by iteratively adjusting key parameters, which can ensure the feasibility of the solution.

The remainder of this paper is organized as follows: Section 2 defines the RE-PFD, while Section 3 proposes the operational constraints of RE-PFD and establishes a TRSEP model. Section 4 designs a customized linearization correction strategy based on the McCormick method. In Section 5, the numerical results of test systems are analyzed. Finally, we draw conclusions in Section 6.

In this section, we define the RE-PFD by referring to the definition of CEF [12]. CEF considers carbon emissions as a parameter in the power flow, thereby achieving the tracking of carbon emissions from all the components. Similarly, RE power can be regarded as a parameter in the power flow.

For convenience, we omit the subscripts scenario s and time t. Taking transmission line l as an example, we define $P_l^{\mathrm{L},\mathrm{RE}}$ as the RE power that passes through line l. Then, the ratio of $P_l^{\mathrm{L},\mathrm{RE}}$ to the line power $P_l^{\mathrm{L}}$ can be defined as RE-PFD, which can be expressed as ${\rho }_l^{\mathrm{L},\mathrm{RE}}$ in Eq. (1).

Suppose that the units located at bus b are regarded as lines with positive power injection; then, the RE-PFD of bus b can be expressed as follows:

where ${\Omega }_{\mathrm{L}\left(b\right)}^{+}$ is the set of lines with a positive power injection connected to bus b.

According to Eqs. (1), (2), the meaning of RE-PFD can be clearly understood, where Eq. (1) explains that the proportion of power provided by RE units in the total power passing through line l is ${\rho }_l^{\mathrm{L},\mathrm{RE}}$, and Eq. (2) explains that the proportion of power provided by RE units in the total injected power of bus b is ${\rho }_b^{\mathrm{B},\mathrm{RE}}$.

For $m\in {\Omega }_{\mathrm{L}\left(b\right)}^{-}$ (where ${\Omega }_{\mathrm{L}\left(b\right)}^{-}$ is the set of lines with a positive power outflow connected to bus b) and $n\in {\Omega }_{\mathrm{L}\left(b\right)}^{+}$, as shown in Fig. 1, according to the principle of proportional sharing [11], the power component of line n passing through line m, denoted as $P_{m\left(n\right)}^{\mathrm{L}}$, satisfies the following relationship:

The RE power that passes through line m can be expressed as Eq. (4).

Combining Eqs. (1), (2), (3), (4), we can deduce the form of Eq. (5), which represents the RE-PFD of line m. It can be observed that the RE-PFD is similar to the CEF density and meets the requirement that the RE-PFD of line $m\in {\Omega }_{\mathrm{L}\left(b\right)}^{-}$ is equal to the RE-PFD of bus b.

Furthermore, according to the definition, we can determine that the RE-PFD of conventional unit g meets 0 and the RE-PFD of RE unit r meets 1. For ESS unit e, the RE-PFD is a variable that changes with time and must be subject to relevant operational constraints. A small example of RE-PFD is shown in Fig. 2.

In this section, based on the definition of RE-PFD, we propose the operational constraints of RE-PFD and establish a TRSEP model for renewable-dominant power systems to cope with the impact of climate change and achieve the development target of RE penetration.

Since the RE-PFD of line l is related to the actual power flow direction, the direction variable ${\tau }_{l,s,t}^{\mathrm{L}}$ is introduced, and the power flow is represented by two non-negative variables $P_{l,s,t}^{\mathrm{L},+}$ and $P_{l,s,t}^{\mathrm{L},-}$. The relationships between this variables can be expressed by Eqs. (6), (7), (8), where $P_{l,s,t}^{\mathrm{L}}$ is operation power of line l at time t in scenario s; $P_{l,s,t}^{\mathrm{L},+}$ and $P_{l,s,t}^{\mathrm{L},-}$ are forward and backward operation power of line l at time t in scenario s, respectively. ${\Omega }_{\mathrm{L}}$ and ${\Omega }_{\mathrm{S}}$: the sets of all transmission and representative scenario, respectively. ${\Omega }_{\mathrm{T}}$: the set of time periods in each scenario, from 1 to |T|; in this paper, a typical day represents a scenario, and |T| is set to 24. ${\tau }_{l,s,t}^{\mathrm{L}}$: the binary direction variable of line l at time t in scenario s: 1 if forward and 0 otherwise. ${\overline{P}}_l^{\mathrm{L}}$: the maximum operation power of line l. ${\rho }_{l,s,t}^{\mathrm{L},\mathrm{RE}}$: the RE-PFD of line l at time t in scenario s. Eqs. (7), (8) restrict the RE-PDF of line l, where p(l) and q(l) are the indices of the start and end buses of line l. When the positive power flows from p(l) to q(l), ${\tau }_{l,s,t}^{\mathrm{L}}$ is equal to 1; otherwise, ${\tau }_{l,s,t}^{\mathrm{L}}$ is equal to 0.

The RE-PFD of bus b can be expressed by Eq. (11), which is calculated based on all the injection power of bus b. ${\Omega }_{\mathrm{L},\mathrm{f}\left(b\right)}$ and ${\Omega }_{\mathrm{L},\mathrm{t}\left(b\right)}$ are the sets of lines starting from and ending at bus b, while ${\Omega }_{\mathrm{G}\left(b\right)}$, ${\Omega }_{\mathrm{RE}\left(b\right)}$, and ${\Omega }_{\mathrm{ES}\left(b\right)}$ are the sets of conventional units, RE units, and ESSs, respectively, located at bus b.

where ${\rho }_{b,s,t}^{\mathrm{B},\mathrm{RE}}$ is RE-PFD of bus b at time t in scenario s; $P_{g,s,t}^{\mathrm{G}}$ and $P_{r,s,t}^{\mathrm{RE}}$ are the operation power of conventional unit g and RE unit r at time t in scenario s, respectively; ${\rho }_{e,s,t}^{\mathrm{ES},\mathrm{RE}}$ is the RE-PFD of ESS unit e at the end of time t in scenario s; $P_{e,s,t}^{\mathrm{ES},\mathrm{d}}$ is discharging power of ESS unit e at time t in scenario s; and ${\Omega }_{\mathrm{B}}$ is the set of buses.

The RE-PFD of ESS unit e is related to the stored RE energy; thus, the form of the RE-PFD constraints is similar to that of the operational constraints [17]. Eq. (12) describes the relationship among the stored RE energy, operation power, and RE-PFD. Eq. (13) represents the lost RE energy of ESSs, and Eq. (14) is the expression of RE-PFD. Eq. (15) limits the initial stored RE energy.

where b(e) is the index of the bus connected to ESS unit e. $E_{e,s,t}^{\mathrm{ES},\mathrm{RE}}$ and $E_{e,s,t}^{\mathrm{ES}}$ are the stored RE energy and total energy of ESS unit e at the end of time t in scenario s, respectively. It should be noted that $E_{e,s,t}^{\mathrm{ES},\mathrm{RE}}$ and ${\rho }_{e,s,t}^{\mathrm{ES},\mathrm{RE}}$ are the instantaneous values at the end of time t. $P_{e,s,t}^{\mathrm{ES},\mathrm{c}}$ is charging power of ESS unit e at time t in scenario s; ${\eta }_e^{\mathrm{ES},\mathrm{c}}$ and ${\eta }_e^{\mathrm{ES},\mathrm{d}}$ are charging and discharging efficiency of ESS unit e, respectively. $E_{e,s,0}^{\mathrm{ES}}$ and $E_{e,s,0}^{\mathrm{ES},\mathrm{RE}}$ are initial stored total energy and RE energy of ESS unit e in scenario s, respectively; ${\rho }_{e,s,0}^{\mathrm{ES},\mathrm{RE}}$ is the initial RE-PFD of ESS unit e in scenario s. ${\Omega }_{\mathrm{ES}}$ is the set of ESSs. Therefore, when ESS unit e is discharging at time t in scenario s, the RE energy released during time t is ${\rho }_{e,s,t-1}^{\mathrm{ES},\mathrm{RE}}{P}_{e,s,t}^{\mathrm{ES},\mathrm{d}}/{\eta }_e^{\mathrm{ES},\mathrm{d}}$.

Eqs. (14), (15), (16) are important constraints related to the key parameters ${\beta }^{\mathrm{R}\mathrm{E}}$, ${\gamma }_o^{\mathrm{O},\mathrm{R}\mathrm{E}}$, and ${\beta }_d^{\mathrm{D},\mathrm{R}\mathrm{E}}$. It can be observed that the actual power generation of RE includes three parts, two of which are transmitted to internal loads and external loads, while the remaining part is lost by ESSs. Eq. (16) specifies that the RE energy supplied to the internal loads should meet the requirement, where ${\beta }^{\mathrm{R}\mathrm{E}}$ represents the minimum ratio of the RE energy to system internal load energy. In Eq. (16), the RE energy consumed by the total internal loads can be obtained by subtracting the RE energy consumed by external loads and the lost RE energy of ESSs from the total RE generation. Eq. (17) limits the proportion of RE energy to meet the requirements for each external load, where ${\gamma }_o^{\mathrm{O},\mathrm{R}\mathrm{E}}$ represents the minimum ratio of RE energy to the energy of the system external load o. Eq. (18) limits the proportion of RE energy to meet the requirements for each internal load, where ${\beta }_d^{\mathrm{D},\mathrm{R}\mathrm{E}}$ represents the minimum ratio of RE energy to the energy of internal load d. It can be observed that parameters ${\gamma }_o^{\mathrm{O},\mathrm{R}\mathrm{E}}$ and ${\beta }_d^{\mathrm{D},\mathrm{R}\mathrm{E}}$ will affect the changes in the transmission path and the geographical distribution of RE. $p_s$ is the weight of scenario s; ${\rho }_{b\left(o\right),s,t}^{\mathrm{B},\mathrm{RE}}$ and ${\rho }_{b\left(d\right),s,t}^{\mathrm{B},\mathrm{RE}}$ are the RE-PFD of the buses where the external load o and the internal load d are located; $P_{d,s,t}^{\mathrm{D}}$ and $P_{o,s,t}^{\mathrm{O}}$ are demand power of system internal load d and external load o at time t in scenario s; ${\Omega }_{\mathrm{D}}$ and ${\Omega }_{\mathrm{O}}$ are the sets of system internal loads and external loads, respectively.

Since Eqs. (9), (10), (11), (12) contain many bilinear variables, it is necessary to find an appropriate linearization method. Furthermore, the focus of this model is to calculate the transmission path of RE, which can appropriately simplify other operational constraints to reduce the complexity of the model. Therefore, this paper adopts direct current (DC) power flow and the economic dispatch model in the operation stage to reduce the computational complexity. On this basis, the network constraints can easily be replaced with the form of alternating current (AC) power flow, and the operation stage can be extended to the unit commitment model. These factors will not affect the calculation method of the RE transmission path.

The detailed form of the TRSEP model is as follows; it contains an objective function, investment constraints, and operational constraints [22]. Eqs. (6), (7), (8), (9), (10), (11), (12), (13), (14), (15), (16) should also be considered in the operational constraints. In addition, typical scenarios are adopted to describe the uncertainties of RE and load power affected by climate change.

For the construction of ESSs, in order to make the model of the ESSs more universal, this section models ESSs uniformly. The set of ESSs can include various types of ESSs, and the differences in their specific operating mechanisms are mainly reflected in the setting of different parameters, such as the ratio of power and energy, the upper and lower limits of stored energy, the regulation period, and the charging/discharging efficiency. Therefore, in practical planning work, planning decision-makers can use a unified mathematical model of ESSs to plan for different types of ESSs by setting different parameters.

The objective function, Eq. (19), minimizes the sum of the annual investment costs and the expected annual operation costs. Eqs. (18), (19), (20) express the annual investment costs of the RE units IC^RE, transmission lines IC^L, and ESSs IC^ES, respectively. Eq. (23) expresses the expected annual operation costs OC.

where $c_r^{\mathrm{RE},\mathrm{I}}$ is the annual power investment cost of RE unit r; $n_r^{\mathrm{RE}}$ is the expanded power of RE unit r; $n_l^{\mathrm{L}}$ is the binary decision variable for the construction of line l: 1 if constructed and 0 otherwise. $c_l^{\mathrm{L}}$ is the annual investment cost of transmission line l; $c_e^{\mathrm{ES},\mathrm{E}}$ and $c_e^{\mathrm{ES},\mathrm{P}}$ are annual energy and power investment cost of ESS unit e, respectively. $n_e^{\mathrm{ES},\mathrm{P}}$ and $n_e^{\mathrm{ES},\mathrm{E}}$ are expanded power and energy of ESS unit e, respectively. K is conversion factor of scenario operation costs; in this paper, K is set to 365 to convert scenario operation costs to expected annual costs. $c_g^{\mathrm{G}}$ is the operation cost of conventional unit g; and $P_{g,s,t}^{\mathrm{G}}$ is operation power of conventional unit g at time t in scenario s.

Eq. (24) is the construction status constraint of the existing lines. Eqs. (23), (24) are the expanded power constraints of the RE units and ESSs. Eq. (27) restricts the relationship between the expanded power and energy of the ESSs.

where ${\Omega }_{\mathrm{L}}^{\mathrm{c}}$ is the set of all candidate lines; ${\overline{n}}_r^{\mathrm{RE}}$ and ${\underline{n}}_r^{\mathrm{RE}}$ are maximum allowable and existing installed capacity of RE unit r, respectively; ${\overline{n}}_e^{\mathrm{ES},\mathrm{P}}$ and ${\underline{n}}_e^{\mathrm{ES},\mathrm{P}}$ are maximum allowable and existing installed capacity of ESS unit e, respectively; and ${\gamma }_e^{\mathrm{ES}}$ is the ratio of expanded energy to expanded power of ESS unit e.

Eqs. (26), (27) are the output range constraint and ramp rate constraint of conventional units, respectively, and Eq. (30) is the output range constraint of the RE units. Eqs. (29), (30), (31), (32), (33), (34) are the operational constraints of the ESSs. Eq. (31) is the mutually exclusive constraint of operational statuses. Eqs. (30), (31) limit the range of the charging and discharging power, respectively. Eq. (34) limits the range of stored energy. Eq. (35) restricts the stored energy at the beginning and ending of scenario s. Eq. (36) describes the intertemporal relationship among the charging power, discharging power, and stored energy.

where ${\overline{P}}_g^{\mathrm{G}}$ and ${\underline{P}}_g^{\mathrm{G}}$ are maximum and minimum operation power of conventional unit g, respectively; ${\Omega }_{\mathrm{G}}$ is the set of conventional units; $\text{R}{\text{U}}_g$ and $\text{R}{\text{D}}_g$ are maximum ramp-up and ramp-down rate of conventional unit g, respectively; ${\alpha }_{r,s,t}^{\mathrm{RE}}$ is the ratio of resource power of RE unit r at time t in scenario s; $x_{e,s,t}^{\mathrm{ES},\mathrm{c}}$ and $x_{e,s,t}^{\mathrm{ES},\mathrm{d}}$ are binary charging and discharging status variables of ESS unit e at time t in scenario s, respectively; ${\underline{\alpha }}_e^{\mathrm{ES}}$ is the minimum stored energy ratio of ESS unit e; ${\underline{n}}_e^{\mathrm{ES},\mathrm{E}}$ is the existing installed energy of ESS unit e; and ${\alpha }_{e,s,0}^{\mathrm{ES}}$ is the initial stored energy ratio of ESS e in scenario s.

Eq. (37) is the power balance constraint of each bus, where ${\Omega }_{\mathrm{D}\left(b\right)}$ and ${\Omega }_{\mathrm{O}\left(b\right)}$ are the sets of the system’s internal loads and external loads located at bus b. Eq. (38) limits the range of the phase angle of the buses. Eq. (39) limits the phase angle of the reference bus b_ref to 0 at any time. Eqs. (38), (39) are the DC power flow constraint and capability constraint of the transmission lines, respectively.

where ${\overline{\theta }}_b$ and ${\underline{\theta }}_b$ are maximum and minimum phase angle of bus b, respectively; ${\theta }_{b,s,t}$ is the phase angle of bus b at time t in scenario s; ${\theta }_{p\left(l\right),s,t}$ and ${\theta }_{q\left(l\right),s,t}$ are the phase angles of the starting and ending buses of line l at time t in scenario s, respectively; $x_l$ is the reactance of transmission line l; and ${\overline{P}}_l^{\mathrm{L}}$ is the maximum operation power of line l.

The TRSEP model is a mixed-integer nonlinear programming model and must be linearized before commercial solvers are applied. In this section, the McCormick method is adopted to linearize the bilinear terms, and a customized feasibility correction strategy is designed to obtain a good feasible solution.

The McCormick method relaxes the bilinear term ah by generating a set of McCormick estimators (MEs) [19], as shown in Fig. 3. The constraints of the MEs are shown in Eqs. (40), (41), (42), (43). It can be observed that no binary variables are generated, so the computational burden does not significantly increase.

where w is the product of a and h; $\overline{a}$ and $\overline{h}$ are the upper bounds of a and h, respectively; $\underline{a}$ and $\underline{h}$ are the lower bounds of a and h, respectively.

Although the MEs do not add binary variables, they still expand the feasible region of the original problem. For example, there is a solution $(a_0,h_0,w_0)\in \left\{(a,h,w)\right|\text{Eqs}.\left(42\right)-\left(45\right)\}$, but $w_0\ne a_0{h}_0$. It can be observed that $(a_0,h_0,w_0)$ is the feasible solution to the linearization problem with MEs, but not the feasible solution to the original problem.

Based on the above analysis, we can make the following proposition.

Proposition 1 Suppose that the feasible region and optimal solution of the original problem $P_0$ are $W_0$ and $({\widehat{\mathit{x}}}_0,{\widehat{\mathit{y}}}_0,{\widehat{\mathit{z}}}_0)\in W_0$, and the feasible region and optimal solution of the linearization problem $P_0$ are $W_0$ and $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{\prime })\in W_0^{\prime }$. Then, we can determine that $W_0\subseteq W_0^{\prime }$.

For the TRSEP model, $\mathit{x}=\{n_r^{\mathrm{R}\mathrm{E}},n_l^{\mathrm{L}},n_e^{\mathrm{E}\mathrm{S},\mathrm{P}},n_e^{\mathrm{E}\mathrm{S},\mathrm{E}}\}$, $\mathit{y}=\{P_{g,s,t}^{\mathrm{G}},P_{r,s,t}^{\mathrm{R}\mathrm{E}},P_{e,s,t}^{\mathrm{E}\mathrm{S},\mathrm{c}},P_{e,s,t}^{\mathrm{E}\mathrm{S},\mathrm{d}},x_{e,s,t}^{\mathrm{E}\mathrm{S},\mathrm{c}},x_{e,s,t}^{\mathrm{E}\mathrm{S},\mathrm{d}},E_{e,s,t}^{\mathrm{E}\mathrm{S}},\mathrm{\Delta }{E}_{e,s,t}^{\mathrm{E}\mathrm{S}},{\theta }_{b,s,t},P_{l,s,t}^{\mathrm{L}},{\tau }_{l,s,t}^{\mathrm{L}}\}$, and $\mathit{z}=\{{\rho }_{b,s,t}^{\mathrm{B},\mathrm{R}\mathrm{E}},{\rho }_{e,s,t}^{\mathrm{E}\mathrm{S},\mathrm{R}\mathrm{E}},{\rho }_{l,s,t}^{\mathrm{L},\mathrm{R}\mathrm{E}},E_{e,s,t}^{\mathrm{E}\mathrm{S},\mathrm{R}\mathrm{E}},\mathrm{\Delta }{E}_{e,s,t}^{\mathrm{E}\mathrm{S},\mathrm{R}\mathrm{E}}\}$. It can be observed that x is the vector of the investment variables, y is the vector of the operational variables related to economic dispatching, and z is the vector of the operational variables related to RE-PFD.

Since we can only obtain $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{\prime })$ by solving $P_0^{\prime }$ and perhaps $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{\prime })\notin W_0$, a feasibility correction strategy must be designed to obtain a feasible solution of $P_0$.

Since bilinear terms are generated by multiplying the variables in y and z, errors will occur in Eqs. (9), (10), (11), (12) after linearization, which will affect Eqs. (14), (15), (16). According to the definition of RE-PFD, if y is known, z can be directly calculated according to the power flow results. Therefore, when we obtain $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{\prime })$ by solving $P_0^{\prime }$, the actual value of $\mathit{z}={\widehat{\mathit{z}}}_0^{″}$ can be calculated using Eqs. (7), (8), (9), (10), (11), (12), (13) on the premise that $\mathit{y}={\widehat{\mathit{y}}}_0^{\prime }$, but ${\widehat{\mathit{z}}}_0^{″}$ may violate Eqs. (14), (15), (16).

Furthermore, the constraints of $P_0$ related to z are Eqs. (7), (8), (9), (10), (11), (12), (13), (14), and the other constraints of $P_0$ are all considered in $P_0^{\prime }$ in the linear form. Therefore, when Eqs. (7), (8), (9), (10), (11), (12), (13), (14) are ignored in $P_0$, $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime })$ is a feasible solution for $P_0$.

Based on the above analysis, if we can ensure that ${\widehat{\mathit{z}}}^{″}\in \left\{\mathit{z}\right|\text{Eqs.}\left(9\right)-\left(15\right),\mathit{y}={\widehat{\mathit{y}}}_0^{\prime }\}$ meets Eqs. (14), (15), (16), then $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{″})$ is a feasible solution of $P_0$.

Next, we discuss the specific strategy. For convenience, we take the parameter ${\beta }^{\mathrm{RE}}$ in Eq. (16) as an example to illustrate.

Proposition 2 Suppose that the value of ${\beta }^{\mathrm{RE}}$ and the feasible region of the linearization problem $P_0^{\prime }$ are ${\beta }_0^{\mathrm{RE}}$ and $W_0^{\prime }$, respectively; the value of ${\beta }^{\mathrm{RE}}$ and the feasible region of the linearization problem $P_1^{\prime }$ are ${\beta }_1^{\mathrm{RE}}$ and $W_1^{\prime }$, respectively; and ${\beta }_1^{\mathrm{RE}}>{\beta }_0^{\mathrm{RE}}$. Then, we can conclude that $W_1^{\prime }\subseteq W_0^{\prime }$.

Proof: Suppose that $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{\prime })\in W_0^{\prime }$ is a feasible solution of $P_0^{\prime }$, and the proportion of RE energy corresponding to $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{\prime })$ is ${\stackrel{\sim }{\beta }}_0^{\mathrm{RE}}$; thus, ${\stackrel{\sim }{\beta }}_0^{\mathrm{RE}}>{\beta }_0^{\mathrm{RE}}$. Similarly, suppose that $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{\prime })\in W_1^{\prime }$ is a feasible solution of $P_1^{\prime }$, and the proportion of RE energy corresponding to $({\widehat{\mathit{x}}}_0^{\prime },{\widehat{\mathit{y}}}_0^{\prime },{\widehat{\mathit{z}}}_0^{\prime })$ is ${\stackrel{\sim }{\beta }}_1^{\mathrm{RE}}$; thus, ${\stackrel{\sim }{\beta }}_1^{\mathrm{RE}}>{\beta }_1^{\mathrm{RE}}$. Obviously, we can determine that ${\stackrel{\sim }{\beta }}_1^{\mathrm{RE}}>{\beta }_1^{\mathrm{RE}}>{\beta }_0^{\mathrm{RE}}$. Thus, $({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime },{\widehat{\mathit{z}}}_1^{\prime })$ is also a feasible solution for $P_0^{\prime }$; that is, for $\forall ({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime },{\widehat{\mathit{z}}}_1^{\prime })\in W_1^{\prime }$, $({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime },{\widehat{\mathit{z}}}_1^{\prime })\in W_0^{\prime }$ is satisfied. Therefore, $W_1^{\prime }\subseteq W_0^{\prime }$.

The parameters ${\gamma }_o^{\mathrm{O},\mathrm{R}\mathrm{E}},\forall o$ and ${\beta }_d^{\mathrm{D},\mathrm{R}\mathrm{E}},\forall d$ can be proven in the same way. Furthermore, we can easily extend this argument to Proposition 3.

Proposition 3. Suppose that the optimal objective values of the problem $P_0^{\prime }$ and $P_1^{\prime }$ are ${\hat{\varphi }}_0^{\text{'}}$ and ${\hat{\varphi }}_1^{\text{'}}$, respectively. Since $W_1^{\prime }\subseteq W_0^{\prime }$, for the minimization problems, we can conclude that ${\varphi }_1^{\prime }\ge {\varphi }_0^{\prime }$, as shown in Fig 4.

Therefore, the optimal objective value of the linearization problem will increase as the parameters rise. Suppose that $({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime },{\widehat{\mathit{z}}}_1^{\prime })\in W_1^{\prime }$ is the optimal solution to the linearization problem $P_1^{\prime }$ and ${\widehat{\mathit{z}}}_1^{″}$ is the actual value of z obtained by recalculating Eqs. (7), (8), (9), (10), (11), (12), (13) on the premise that $(\mathit{x},\mathit{y})=({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime })$. According to Proposition 2, we can increase the value of ${\beta }^{\mathrm{RE}}$ to narrow the feasible region $W_1^{\prime }$ of the problem $P_1^{\prime }$. Since $W_1^{\prime }\subseteq W_0^{\prime }$ and $W_0\subseteq W_0^{\prime }$, when ${\beta }^{\mathrm{RE}}$ is adjusted to the appropriate values, $W_1^{\prime }$ can be narrowed to make $({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime },{\widehat{\mathit{z}}}_1^{″})\in W_0$ hold, as shown in Fig. 5. Therefore, $({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime },{\widehat{\mathit{z}}}_1^{″})$ is the feasible solution for $P_0$.

After solving the linearization problem $P_1^{\prime }$ and obtaining the optimal solution $({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime })$, we construct the following subproblem shown in Eqs. (45), (46), (47), (48), (49), (50), (51), (52), (53), (54), (55), (56) to obtain ${\widehat{\mathit{z}}}_1^{″}$ and verify whether ${\widehat{\mathit{z}}}_1^{″}$ is a feasible solution for the original problem $P_0$.

where $\psi$ is the objective function value; $\mathrm{\Delta }{\beta }^{\mathrm{RE}}$, $\mathrm{\Delta }{\gamma }_o^{\mathrm{O},\mathrm{RE}}$, and $\mathrm{\Delta }{\beta }_d^{\mathrm{D},\mathrm{RE}}$ are the slack variables of ${\beta }^{\mathrm{RE}}$, ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$, and ${\beta }_d^{\mathrm{D},\mathrm{RE}}$, respectively; and all variables with the superscript “^” are constants. This subproblem is pure linear programming, which is very easy to solve. If the optimal value $\widehat{\psi }=0$, then ${\widehat{\mathit{z}}}_1^{″}$ satisfies the Eqs. (14), (15), (16) of $P_0$, and $({\widehat{\mathit{x}}}_1^{\prime },{\widehat{\mathit{y}}}_1^{\prime },{\widehat{\mathit{z}}}_1^{″})$ is the feasible solution for $P_0$; otherwise, ${\widehat{\mathit{z}}}_1^{″}$ is infeasible for $P_0$, and it is necessary to continue to adjust the parameters.

The calculation process to solve the TRSEP model problem with the customized feasibility correction strategy is shown in Algorithm 1. For the sake of simplicity, Table 1 illustrates the calculation process using parameters ${\beta }^{\mathrm{RE}}$ and ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ as examples, and parameter ${\beta }_d^{\mathrm{D},\mathrm{RE}}$ can also be iterated according to this process. The process can be divided into two procedures. Procedure 1 (P1) is used to iteratively increase the parameters to obtain a feasible solution, while procedure 2 (P2) adopts the dichotomy to adjust the parameters in order to obtain a better feasible solution. Since the optimal value ${\widehat{\varphi }}_{\left(w\right)}$ is only related to ${\widehat{\mathit{x}}}_{\left(w\right)}$ and ${\widehat{\mathit{y}}}_{\left(w\right)}$, ${\widehat{\varphi }}_{\left(w\right)}$ will not be affected when ${\widehat{\mathit{z}}}_{\left(w\right)}$ is modified to ${\widehat{\mathit{z}}}_{\left(w\right)}^{\prime }$. According to Proposition 3, during the iteration process, we can update the upper and lower bounds by adjusting the parameters.

It can be observed that, in Algorithm 1, when $({\widehat{\mathit{x}}}_{\left(w\right)},{\widehat{\mathit{y}}}_{\left(w\right)},{\widehat{\mathit{z}}}_{\left(w\right)}^{\prime })$ is feasible for $P_0$, ${\widehat{\varphi }}_{\left(w\right)}$ can provide an upper bound; otherwise, ${\widehat{\varphi }}_{\left(w\right)}$ can provide a lower bound. P1 can obtain the first feasible solution, and P2 can reduce the gap between the upper bound (UB) and the lower bound (LB) iteratively by adjusting the parameters. Therefore, a better feasible solution can be obtained by P2.

All calculations are performed using CPLEX 12.80 in the C++ environment on a workstation computer equipped with an Intel(R) Xeon(R) Gold 6136 3.00 GHz processor and 192 GB of random access memory (RAM).

In this section, we select Garver’s six-bus test system and the Institute of Electrical and Electronics Engineers (IEEE)-118 test system to conduct numerical analyses. For all the following cases, we set ${\gamma }_e^{\mathrm{ES}}=2$ and ${\eta }_e^{\mathrm{ES},\mathrm{c}}={\eta }_e^{\mathrm{ES},\mathrm{d}}=95\%$. To facilitate the comparison of the impact of the key parameters ${\beta }^{\mathrm{RE}}$ and ${\gamma }_o^{\mathrm{O},\mathrm{RE}},\forall o$, only wind units are allowed to be constructed, and the investment cost is 550 USD·kW⁻¹ with a 30-year lifetime. The transmission investment cost is 6200 USD·(km·MW)⁻¹ with a 40-year lifetime, and the power and energy investment costs of the ESSs are 211 and 189 USD·(kW·h)⁻¹, respectively, with a 10-year lifetime [23]. The convergence gap is set to 0.1%, and all investment costs are converted into annual investment costs in the target year by Eq. (60), where $\stackrel{\sim }{c}$ represents the annual investment cost, c₀ represents the whole-cycle investment cost, Y represents the lifetime, and a = 0.08 represents the interest rate.

Garver’s six-bus test system contains six buses and six existing transmission lines. We allow all the corridors to expand their lines, with the maximum number of lines in each corridor being three. B2, B4, B5, and B6 are selected to construct ESSs, and the maximum allowable capacity is 100 MW for each bus. All the buses are allowed to construct RE units, and the maximum allowable capacity is 400 MW for each bus. The existing network topology, power source installations, and load levels are shown in Fig. 6. Detailed data for this system is provided in Ref. [24].

In this subsection, we set B2 (Case 1), B4 (Case 2), B5 (Case 3), and B6 (Case 4) as the external load buses (ELBs). The external load capacity is 150 MW. The total internal load energy and external load energy are 4.2664 × 10⁶ and 8.103 × 10⁵ MW·h, respectively, for all the cases. ${\beta }^{\mathrm{RE}}$ is set to 40% and ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ is set to 50%. The cost results, investment results, and RE transmission path results are respectively shown in Table 1, Table 2, Table 3. Here, TC represents the total costs, E^RE is the total RE energy, E^D,RE is the RE energy supplying the internal loads, E^O,RE is the RE energy supplying the external loads, ΔE^ES,RE is the lost RE energy of the ESSs, ${\widehat{\beta }}^{\mathrm{RE}}$ is the actual ratio of E^D,RE to the internal load energy, and ${\widehat{\gamma }}_o^{\mathrm{O},\mathrm{RE}}$ is the actual ratio of E^O,RE to the external load energy.

It can be observed from Table 2 that—driven by the RE development goals—the total installed capacity of RE in all cases exceeds that of the conventional units, reaching over 800 MW. In this situation, ESSs exceeding 150 MW are constructed to ensure power system balance. After testing, it was found that, without considering ESS construction, relying solely on the expansion of transmission lines and RE cannot ensure power system balance, resulting in serious load-shedding issues. Therefore, to cope with the impact of climate change on renewable-dominant power systems, the construction of ESSs is essential.

When the location of the ELB differs, the difficulty of ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ meeting the requirement varies, resulting in differences in the location and capacity of the expanded lines, RE units, and ESSs. Therefore, the cost results will also vary. For example, the conventional capacity near B5 is high, whereas the conventional capacity near B6 is low. Therefore, the requirement of ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ is easier to meet for B6 than for B5. From a cost perspective, the total cost of Case 4 is lower than that of Case 3. As shown in Table 3, when ${\widehat{\gamma }}_o^{\mathrm{O},\mathrm{RE}}$ reaches 50% of B5, ${\widehat{\beta }}^{\mathrm{RE}}$ has exceeded 40%, and when ${\widehat{\beta }}^{\mathrm{RE}}$ reaches 40% of B6, ${\widehat{\gamma }}_o^{\mathrm{O},\mathrm{RE}}$ has exceeded 50%.

For Case 1, the total load capacity of B2 is 390 MW. In addition to the 200 MW of RE capacity at B6, 305 MW of RE capacity is expanded at B2 to meet the requirement of ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$. Meanwhile, the conventional unit at B3 also needs to transmit power to B2 to ensure a power balance. Therefore, B2 and B3 and B2–B6 are constructed. To ensure the internal load demands, RE units are also constructed at B4 and B5. With the connection between B2 and other buses, there are sufficient flexible resources available; thus, less ESS capacity is constructed at B2. Fig. 7 shows the results of the RE transmission path for Case 1.

For Case 2, the total load capacity of B4 is 310 MW. To meet the requirement of ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$, 380 MW of RE capacity is installed at B4, and 100 MW·(200 MW·h)⁻¹ of ESS capacity is installed to mitigate the fluctuation in RE power. Since the large capacity of RE units can supply the loads at B2 and B4, the conventional units at B1 and B3 can mainly transmit power to B5, resulting in a smaller RE capacity at B5. Compared with Case 1, the ESS capacity at B5 is also reduced.

For Case 3, the total load capacity of B5 is 390 MW. As the buses connected to B5 are equipped with conventional units, the RE capacity at B5 has reached the maximum value to meet the requirement of ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$. Meanwhile, it is still necessary to construct 144 MW of RE capacity at B3 to transmit RE power to B5, thereby expanding B3–B5. Correspondingly, B2 and B4 can reduce the expanded capacity of RE and ESS, still meeting the requirement of ${\beta }^{\mathrm{RE}}$.

For Case 4, the total load capacity of B6 is 150 MW. Since the RE capacity is greater than the load capacity at B6, the requirement of ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ is easy to meet. Meanwhile, B3–B6 are constructed so the conventional unit at B3 can cope with the fluctuations in load and RE power at B6, and the ESS capacity at B6 is relatively small. Moreover, high capacity of RE units and ESSs can be constructed at B2 and B5 to meet the internal load requirements.

Based on the above analysis, the decision-making method for renewable-dominated power systems proposed in this paper can effectively distinguish the RE transmission path, thereby meeting both the RE development goals within the system and the requirements for green power trading between systems during the planning stage.

In this subsection, we present the iterative process of the algorithm for Case 1 and Case 2. The iterative parameter values are shown in Table 4, Table 5. The convergence process is shown in Fig. 8, Fig. 9.

As shown in Table 4, when w = 1, both ${\beta }^{\mathrm{RE}}$ and ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ do not meet the requirements; thus, they both need to be increased in P1. When w = 4, the value of ${\widehat{\psi }}_{\left(w\right)}$ is 0, and a feasible solution is obtained. Then, the iterative process will enter P2 at w = 5. In P2, Algorithm 1 adjusts ${\beta }^{\mathrm{RE}}$ and ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ through dichotomy to find a better feasible solution. Here, w = 3 and 4 respectively determine the LB and UB of the parameters in P2. As shown in Fig. 8, LB is updated when an infeasible solution is obtained, and UB is updated when a feasible solution is obtained, where the bold numbers represent the value of ${\widehat{\varphi }}_{\left(w\right)}$ for each iteration. When the gap meets the requirement, the latest feasible solution (w = 8) will be output.

For Case 2, the conventional capacity near B4 is relatively small, and the requirement of ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ is easier to meet than that of ${\beta }^{\mathrm{RE}}$. As shown in Table 5, ${\gamma }_o^{\mathrm{O},\mathrm{RE}}$ remains unchanged, whereas ${\beta }^{\mathrm{RE}}$ needs to increase in P1. When w = 2, a feasible solution is obtained, and the iterative process will enter P2 at w = 3. When the gap meets the requirement, the latest feasible solution (w = 5) will be output.

Moreover, the number of iterations for Cases 3 and 4 is 11 and 5, respectively. It can be observed that the number of iterations is greater than 1 for each case, where w = 1 represents the McCormick method without the feasibility correction strategy. This result indicates that, without the feasibility correction strategy, the feasibility of the solution cannot be ensured at w = 1, demonstrating the necessity of the feasible correction strategy.

On the basis of the parameter settings in Case 4 (where B6 is the ELB), this section further expands the analysis by adding energy transition areas within the system. The parameter ${\beta }_d^{\mathrm{D},\mathrm{RE}}$ is set to 55%. The comparative cases include the following:

Case 5: Set B4 as the energy transition area;

Case 6: Set B5 as the energy transition area;

Case 7: Simultaneously set B4 and B5 as the energy transition areas.

The calculation results for RE in the above cases are shown in Table 6, Table 7.

It can be observed that, when the energy transition area is located at B4 or B5, the planning scheme for Case 4 cannot meet the requirement of ${\beta }_d^{\mathrm{D},\mathrm{RE}}$, and the relevant constraint needs to be set for further planning.

For Case 5, when B4 is selected as the energy transition area, compared with Case 4, the RE capacity is mainly transferred from B2 to B4 and is used to directly provide RE energy to B4. The RE capacity located at B5 remains basically unchanged, and the corresponding change in the value of ${\beta }_d^{\mathrm{D},\mathrm{RE}}$ at B5 is also relatively small. For Case 6, when B5 is chosen as the energy transition area, the RE capacity is mainly transferred from B2 and B4 to B5, and the value of ${\beta }_d^{\mathrm{D},\mathrm{RE}}$ at B4 decreases significantly. For Case 7, the requirements of ${\beta }_d^{\mathrm{D},\mathrm{RE}}$ at B4 and B5 are considered simultaneously, and the RE capacity at B2 is transferred to B4 and B5, in comparison with Case 4.