Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Frequency-constrained Co-planning of Generation and Energy Storage with High-penetration Renewable Energy  PDF

  • Chengming Zhang
  • Lu Liu
  • Haozhong Cheng
  • Dundun Liu
  • Jianping Zhang
  • Gang Li 2
Key Laboratory of Control of Power Transmission and Conversion, Shanghai Jiao Tong University, Shanghai, China; East China Branch of State Grid Corporation, Shanghai, China

Updated:2021-08-02

DOI:10.35833/MPCE.2020.000743

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

Large-scale renewable energy integration decreases the system inertia and restricts frequency regulation. To maintain the frequency stability, allocating adequate frequency-support sources poses a critical challenge to planners. In this context, we propose a frequency-constrained coordination planning model of thermal units, wind farms, and battery energy storage systems (BESSs) to provide satisfactory frequency supports. Firstly, a modified multi-machine system frequency response (MSFR) model that accounts for the dynamic responses from both synchronous generators and grid-connected inverters is constructed with preset power-headroom. Secondly, the rate-of-change-of-frequency (ROCOF) and frequency response power are deduced to construct frequency constraints. A data-driven piecewise linearization (DDPWL) method based on hyperplane fitting and data classification is applied to linearize the highly nonlinear frequency response power. Thirdly, frequency constraints are inserted into our planning model, while the unit commitment based on the coordinated operation of the thermal-hydro-wind-BESS hybrid system is implemented. At last, the proposed model is applied to the IEEE RTS-79 test system. The results demonstrate the effectiveness of our co-planning model to keep the frequency stability.

I. Introduction

AS serious energy crisis spreads in the world, renewable energy sources (RESs) are deemed as the best way to achieve the sustainable development [

1]. Besides the deployment of large-scale power plants in the transmission network, RESs have also been efficiently utilized in the distribution network, microgrid [2], and even ship system [3] nowadays. However, RES-based power plants are connected to the power grid through converters, which decouple the frequency from the speed of rotating machines. Thus, these electronic devices neither provide inertia response nor participate in the auxiliary service market. With high penetration of renewable energy in the future, the extensive replacement of synchronous generators (SGs) decreases the system inertia and weakens the frequency regulation. Frequency stability has recently attracted the interest of researchers.

Nowadays, massive frequency response modes are excavated to supply adequate system inertia and maintain frequency within the acceptable range. RES-based power plants are studied to obtain virtual inertia and primary frequency response (PFR) by integrating frequency control loops. In this respect, various control strategies such as deloading and overproduction have been developed to facilitate the provision of frequency support from variable speed wind turbines (VSWTs) [

4]-[7]. In addition, the battery energy storage system (BESS) can also respond to the power imbalance within milliseconds through fast-dependent frequency responses [8]. In order to quickly stabilize the frequency of different areas, the tie-line power is adopted as the additional input signal of demand response [9]. However, the impacts of synthetical frequency responses on both operation and planning problems are not yet fully understood.

In recent years, frequency-constrained problems mainly concentrate on optimal operation such as economic dispatch (ED) and unit commitment (UC). For instance, [

10] inserts nonlinear frequency deviation constraints into AC optimal power flow (ACOPF) to keep the frequency stability. Considering credible generation outages, [11] takes full usage of primary frequency reserves to satisfy the minimum frequency deviation. However, unit states are totally neglected in ED problems. In this vein, [12] proposes a mixed-integer linear UC to co-optimize the energy production, inertia-dependent fast frequency responses, and PFR at the same time. Virtual inertia and droop controls are adopted to reduce the frequency support from SGs in [13]. Some researches [14], [15] also design the frequency reserve market to provide ancillary services and maximize benefits. However, dynamic characteristics of PFR from SGs are omitted in such researches, which possibly overstates response speeds and available amounts of PFR. Furthermore, a stochastic UC model that incorporates both SGs and grid-connected inverters to obtain synthetical frequency responses is proposed in [16] to investigate the dynamic frequency in low-inertia power systems. Considering the large-scale wind power participating in frequency responses, [17] proposes a UC model with respect to coordination of steady state and transient state based on an inner-outer iterative optimization method. However, frequency responses from hydro units are rarely considered in existing studies, whereas a fast inertial response and PFR can be timely obtained through hydro turbines and governors [18].

With respect to planning problems, the frequency stability has been included into the siting and sizing issues of BESS. For instance, the contribution of inertia and droop controls from BESS is analyzed in [

19]. An operation-constrained BESS planning model is presented to supply enough primary frequency reserves in [20], [21], whereas the frequency stability cannot be guaranteed due to the basic indices such as frequency nadir, rate-of-change-of-frequency (ROCOF), and quasi steady-state frequency deviation are totally omitted. Furthermore, a frequency-constrained BESS expansion model with respect to frequency nadir is proposed in [22], where a linear ramp speed replaces the dynamic characteristics of primary frequency feedbacks. Also, overlooking frequency responses from RES-based power plants incurs unnecessary BESS allocation. Besides, recent developments in generation expansion planning (GEP) problem mainly highlight the need for RES integration [23], [24], combination with short-term operation [1], system reliability [25], and market design [26]. In the future low-inertia power systems, traditional GEP methods cannot make a trade-off between economy and frequency stability, which possibly causes serious frequency collapse ultimately.

In this paper, we mainly focus on RES planning that installs wind farms to facilitate renewable consumption. RES planners should impose frequency limits to secure frequency stability and avoid system collapse. In practice, few existing studies and projects have considered the frequency stability when making RES planning decisions. Without such considerations, a system planner tends to overinvest in RES, while the total delivery of inertial and PFR may be reduced in real time, leading to serious frequency deviations or renewable production curtailment. Therefore, a frequency-constrained planning model of generation and BESS has been coordinated to simultaneously satisfy the frequency stability and power consumption requirements. In general, the main contributions are threefold.

1) A novel multi-machine system frequency response (MSFR) model that incorporates thermal units, hydro units, wind farms, and BESS is proposed to analyze the dynamic frequency trajectory. On this basis, the multi-machine aggregation method is utilized to aggregate the MSFR model into an equivalent single-generator by converting all control loops to a thermal form. Meanwhile, the available frequency response power and power-headroom constraints are both constructed to activate sufficient PFR.

2) Two frequency constraints, namely the limitation on ROCOF and frequency nadir, are deduced from aggregated MSFR (AMSFR) model. The frequency nadir is further transformed into a system frequency response power constraint (FRPC). Furthermore, a data-driven piecewise linearization (DDPWL) method based on hyperplane fitting and data classification is proposed to linearize the highly nonlinear FRPC by solving a mixed-integer second-order cone problem.

3) A frequency-constrained coordination model of generation expansion planning and energy storage installation (GEP&ESI) is presented to supply enough system inertia and activate PFR. The short-term UC is inserted into a long-term planning procedure considering renewable portfolio standards (RPSs) to obtain high-penetration renewable integration. Meanwhile, a hybrid thermal-hydro-wind-BESS operation model is proposed to promote renewable consumption, which constructs extra water constraints for cascading hydro units. Afterwards, linearization methods comprising the big-M method, reformulation linearization technique (RLT), and McCormick envelopes are adopted to deal with bilinear constraints, which ultimately forms a mixed-integer linear optimization problem. Primary frequency reserves are also confined to satisfy the power-headroom requirement in the proposed MSFR model.

The remainder of this paper is organized as follows. The MSFR model is introduced in Section II. A DDPWL method is proposed and discussed in Section III. Then the GEP&ESI model that incorporates frequency limitations is presented in Section IV. The case studies are conducted and explained in Section V. Section VI concludes this paper.

II. MSFR Model

A. Frequency Responses from Wind Farms

RES-based inverters are always used to supply both inertial and droop controls through virtual synchronous machine (VSM) technique. In brief, the virtual inertia mimics synchronous inertia and releases over-generated power within milliseconds [

6]. However, the maximum power point tracking (MPPT) mode designed for wind farms cannot supply PFR since no power-headroom exists in the face of power imbalance. In this vein, the de-loading operation should be adopted instead, which ensures sufficient spaces to regulate output power based on the rotor speed adjustment.

The frequency response power from wind farms can be expressed as (1) in the time domain. Meanwhile, the time delay is inserted into the transfer function in the frequency domain, as shown in (2). Since the frequency response power depends on de-loading operation, the maximum amounts of PFR are limited by renewable curtailment in constraint (3), which relates to the droop constant and the maximum frequency deviation at the same time.

Δpw(t)=-2HwdΔf(t)dt-DwΔf(t) (1)
Δpw(s)Δf(s)=-2Hws-Dw1+Tws (2)
pw,tcDwpw,tfλwΔfmax (3)

where Δf(t) and Δf(s) are the frequency deviations in the time domain and frequency domain, respectively; Δpw(t) and Δpw(s) are the total frequency response power in the time domain and frequency domain, respectively; pw,tf and pw,tc are the forecasted output power and curtailment power for wind farm w at time t, respectively; Δfmax is the maximum frequency deviation; Hw, Dw, and Tw are the inertia constant, droop constant, and response time of wind farm w, respectively; and λw is the reserve coefficient, which can be used to regulate the preset power-headroom for PFR.

B. Frequency Responses from BESS

In this paper, Lithium-ion batteries are selected as BESS devices, which activate quick PFR within 40 ms and keep the duration time from few minutes to hours [

1]. In general, the dynamic characteristics of exchanging power between BESS and the power grid can be represented by the first-order transfer function as shown in (4).

ΔPe(s)Δf(s)=KeRe(1+Tes) (4)

where ΔPe(s) is the frequency response power of BESS e in the frequency domain; and Ke, Re, and Te are the gain coefficient, droop constant, and response time of BESS e, respectively.

Since the droop control is adopted to supply PFR, the timely frequency response power should be limited within the preset power-headroom. Similarly, the power-headroom depicted in Fig. 1 can guarantee sufficient responses. The left block refers to normal operation, while the right one relates to the frequency regulation. It can be found that BESS charges (point A) or discharges (point B) at the preset power in normal operation. With respect to the positive power imbalance (due to generator contingency or demand increase), BESS is forced to increase the output power for supplying an upward frequency response. This can be obtained through two ways, i.e., ① increasing the discharging power; and ② reducing the charging power even transforming to the discharging state. That means the preset power (points A and B) can be extremely adjusted to point D for obtaining the maximum upward frequency responses. Similarly, the downward frequency response power represent the power-headroom with incremental generation. In Fig. 1, Ee,tc and Ee,td are the upward PFR power of charging and discharging at the beginning of frequency regulation, respectively; Ee,tc' and Ee,td' are the downward PFR power of charging and discharging at the beginning of frequency regulation, respectively; Rec,max and Red,max are the maximum charging and discharging power, respectively; and re,tc and re,td are the realistic charging and discharging power, respectively.

Fig. 1 Schematic diagram of PFR from BESS.

On this basis, constraints (5) and (6) are proposed to stipulate ranges of frequency response power with respect to discharging and charging states. Constraint (7) gives an expression of total PFR power. Only one of Ee,tc and Ee,td would take effect at any particular time, which can be guaranteed by integrating binary indicators (ve,tc and ve,td) and adding exclusive constraint (8). The stored energy after frequency responses should be subject to the permitted energy ranges in (9) and (10), which provides enough energy headroom to sustain frequency responses for the duration time of PFR ΔtP. Note that (9) gives out a clear correlation of the stored energy between the normal operation and post-regulation. Here, parts ①-③ refer to the changeable energy due to the frequency regulation. Specifically, part ① relates to the condition that BESS discharges power in normal operation. BESS charges in both normal operation and frequency regulation with part ②. On the contrary, BESS charges in normal operation but transforms to the discharging state after frequency regulation with part ③. Constraint (10) restricts the exclusive state through binary variables ke,tc and ke,td.

0Ee,td(Red,max-re,td)ve,td (5)
0Ee,tc(Red,max+re,tc)ve,tc (6)
Ee,t=Ee,tc+Ee,td (7)
ve,td+ve,tc1 (8)
SOCeminSOCe,t'=SOCe,t-(Ee,td+re,td)ΔtP/ηed+ηec(re,tc-Ee,tc)ke,tcΔtP-(Ee,tc-re,tc)ke,tdΔtP/ηedSOCemax (9)
ke,tc+ke,td1    ke,tc{0,1},ke,td{0,1} (10)

where ve,td and ve,tc are the indicators of discharging and charging operation states, respectively; Ee,t is the total PFR power; SOCe,t and SOCe,t' are the stored energy before and after frequency responses, respectively; SOCemin and SOCemax are the permitted minimum and maximum stored energy, respectively; and ηed and ηec are the discharging and charging efficiencies of BESS e, respectively.

C. Formulation of AMSFR Model

An innovative MSFR model comprising synthetical frequency responses from SGs (both thermal and hydro) as well as converters (both wind farm and BESS) is integrated to maintain frequency stability. This MSFR model can be illustrated by the block diagram in Fig. 2.

Fig. 2 Schematic diagram of MSFR model.

For thermal units, Fh and Tr are the fraction of power generated by the high-pressure turbine and the reheat time of thermal unit, respectively; and Kg and Rg are the mechanical power gain coefficient and the governor speed constant, respectively. For hydro units, frequency responses mainly come from the synchronous inertia, governors, and hydro turbines, where the core component is the governor [

18]. Thus, turbines are omitted in this paper since a low-hydraulic power system is considered with a negligible water hammer effect. In this vein, RP and RC=RT/RP are the permanent and synthetical droop constants for hydro units, respectively, and RT is the temporary droop constant for hydro units; Td is the governor time of hydro units. Moreover, ΔPG, ΔPH, ΔPW, and ΔPE are the PFR power from the thermal unit, hydro unit, wind farm, and BESS in response to frequency deviation Δf, respectively; ΔPL is the system unbalanced power; and HE and DE are the equivalent system inertia and load damping constant, respectively.

As can be observed from Fig. 2, large-scale control loops for multi-machines and converters incur unacceptable computation obstacles in frequency calculation. Thus, we aggregate massive devices into an equivalent machine by means of the parameter equivalence method from [

27], which can construct an AMSFR model in Fig. 3. Here, FH, TR, and RT are equivalent frequency parameters.

Fig. 3 Schematic diagram of AMSFR model.

Furthermore, the power gain factor Km of each synchronous generator and converter can be defined as:

Km,i=Pi,maxSB,sys    i𝒢T𝒢H𝒳 (11)

where Pi,max is the maximum response capacity of device i; SB,sys is the system base power; and 𝒢T, 𝒢H, , and 𝒳 are the sets of thermal units, hydro units, wind farms, and BESS, respectively.

In this paper, we assume that the system base power equals the hourly demand so that the damping coefficient always keeps the same in any time interval. Thus, the equivalent system inertia can be calculated by (12), which relates to the individual inertia and response capacity at the same time.

HE,t=1bdb,tg𝒢HgPg,maxug,t (12)

where HE,t is the equivalent system inertia at time t; Hg is the inertia constant for generator g; Pg,max is the maximum frequency response capacity of generator g; db,t refers to demand b at time t; ug,t is the frequency-support state of generator g at time t; and 𝒢 and are the sets of SGs and power grid buses, respectively.

To derive the equivalent frequency parameters, the transfer functions of wind farms, BESS, and hydro units can be converted into (13)-(15), respectively, in sequence.

2Hws+Dw1+Tws=Dw1+2HwTws/(TwDw)1+Tws=1Rw1+FwTws1+TwsRw=1DwFw=2HwTwDw (13)
KeRe(1+sTe)=1Re/Ke1(1+sTe)=1Re'1+FeTes1+TesRe'=ReKeFe=0 (14)
1RP1+Tds1+RCTds=1RP1+(1/RC)RCTds1+RCTds=1RP1+FdTd's1+Td'sTd'=RCTdFd=1RC (15)

where Rw and Re' are the equivalent droop coefficients of wind farm and BESS, respectively; Fw, Fe, and Fd are the equivalent power fractions of wind farm, BESS, and hydro unit, respectively; and Td' is the equivalent time constant of hydro unit.

Afterwards, the newly-defined droop and governor speed constants can be described as (16) and (17), respectively. Similarly, the occupation proportion is defined as (18), and all proportions should sum up to one. On this basis, we construct the linear formula of FH and TR in (19) and (20), respectively.

κi=Km,iRi    i𝒢T𝒢H𝒳 (16)
1RT,t=i𝒢T𝒢H𝒳κi=1bdb,tg𝒢TKgRgPg,maxug,t+d𝒢H1RdPd,maxud,t+w1Rwpw,tfuw,t+e𝒳1Re'Remaxue,t (17)
λi=κii𝒢T𝒢H𝒳κii𝒢T𝒢H𝒳λi=1 (18)
FH,tRT,t=1RT,ti𝒢T𝒢HλiFh,i=1bdb,tg𝒢TKgPg,maxFgug,tRg+d𝒢HFdPd,maxud,tRd+wpw,tfFwuw,tRw (19)
TR,t=1bdb,tg𝒢TKgPg,maxTr,gug,tRg+d𝒢HTd'Pd,maxud,tRd+wpw,tfTwuw,tRw+e𝒳KeRemaxTeue,tRe' (20)

where ud,t, uw,t, and ue,t are the frequency-support states of hydro unit, wind farm, and BESS at time t, respectively; Pd,max is the maximum capacity of hydro unit d; and Remax is the maximum response capacity of BESS e.

Based on the block diagram in Fig. 3, the transfer function G(s) of the proposed AMSFR model can be deduced as (21). Here, the natural oscillation frequency wn and damping ratio ξ are accordingly calculated by (22).

G(s)=Δf(s)ΔPL(s)=RTwn2(1+TRs)(DERT+1)(s2+2ξwns+wn2) (21)
wn2=DERT+12RTHETRξ=RTDETR+2RTHE+FHTR2(DERT+1)wn (22)

The analytical formulation of Δf(t) in the time domain can be given as (23), where the damping frequency wr and coefficients α and φ are defined in (24). Then, a step signal is settled for the power deviation.

Δf(t)=L-1G(s)ΔPLs=RTΔPLDERT+1(1+αe-ξwntsin(wrt+φ)) (23)
α=1-2TRξwn+TR2wn21-ξ2wr=wn1-ξ2φ=arctanwrTR1-TRξwn-arctan1-ξ2-ξ (24)

According to the frequency deviation (23), we implement its derivation to deduce the maximum frequency deviation Δfmax and related nadir time tmax. These two metrics expressed as (25) are directly decided by the equivalent frequency parameters of system. In addition, the ROCOF f˙max is employed to limit the maximum changing speed of frequency at the initial time of frequency regulation, which enables to avert unexpected generator tripping. The formulation of ROCOF is shown as (26) [

16], which mainly relates to the system inertia and the magnitude of the unbalanced power.

Δfmax=RTΔPLDERT+11+1-ξ2αe-ξwntmaxtmax=1wrarctanwrTRξwnTR-1 (25)
f˙max=dΔf(t0+)dt=ΔPL2HE (26)

Furthermore, the frequency deviation and ROCOF restrictions are constructed as (27) and (28), respectively. In this vein, the ROCOF limitation equals the system inertia requirement, which is totally linear and can be directly put into the optimization model. Moreover, the frequency deviation is transformed into the power expression in constraint (29) to reflect the basic requirement of system frequency response power. However, the right part defined as PFR(w0) is a highly nonlinear multi-variable function, which is decided by the equivalent frequency parameters w0=[HE,RT,FH,TR,DE]T and should be further linearized in the next section.

fbRTΔPLDRT+11+1-ξ2αe-ξwntmaxΔfmax (27)
HEfbΔPL2ROCOFmax (28)
ΔPL(DRT+1)ΔfmaxRT1+1-ξ2αe-ξwntmax=PFR(w0) (29)

where fb is the nominal frequency; and ROCOFmax is the maximum permitted ROCOF.

III. DDPWL Method

In this section, a traditional piecewise linearization (PWL) method from [

16], [28] is introduced, which is generally applied to linearize the multi-variable functions. In brief, this method maximizes all hyperplanes to approach the accurate function value. Before linearization, two widely-used assumptions are made as follows.

1) Due to the low sensitivity of the frequency deviation to the reheat time constant (about 1%) [

28], the identical value is always assumed for all frequency-support devices without loss of accuracy.

2) With respective damping and droop gains usually strictly prescribed within narrow ranges by the system operator, the equivalent damping constant of the system can be settled as a constant [

16].

On this basis, the decision variable set can be further reduced as win=[HE,RT,FH]T. To form a linear correlation between PFL(w) and frequency-support states, fitting variables coming from (17) and (19) are employed instead. Thus, the improved PWL function can be given as:

PFL(w)=maxn𝒫cnhHT+cnr1RT+cnfFHRT+bn (30)
w=HE,1RT,FHRTT (31)

where cn=[cnh,cnr,cnf] and bn refer to the fitting coefficients, which are optimized to eliminate fitting errors between the real value PFR(w) and the fitting value PFL(w); n counts the number of adopted hyperplanes; and 𝒫 is the set of hyperplanes.

Note that the second-order terms are always optimized to derive the appropriate fitting coefficients based on the dataset 𝒟. The traditional PWL problem can be constructed as:

mincn,bni𝒟(PFL(wi)-PFR(wi))2 (32)

Since the problem (32) minimizes the squared errors in the data set, the frequency stability is possibly destroyed at some points. Meanwhile, the over-conservativeness takes place if a lower fitting value adopts and facilitates more frequency supports. Meanwhile, enough data points and hyperplanes should be inserted to guarantee accuracy so that solution obstacles always exist. To overcome the shortcomings above, [

16] proposes an alternative approach called bound extraction to confine the values within a plausible range, which guarantees that the frequency nadir threshold is not violated. This method introduces linear bound extraction constraints to reduce solution time and guarantee conservativeness for the UC problem. Nevertheless, it should enumerate every possible frequency-support state to take the worst-case frequency parameters as the restricted bounds, where 2NG combinations for all frequency-support conditions nearly lead to an untractable solution in large-scale power systems. Note that NG is the total number of frequency support sources.

In this paper, a modified DDPWL methodology is proposed to deal with the nonlinear FRPC. Firstly, the initial problem (30)-(32) can be transformed into the problem (33)-(35) by defining fL(wi)=maxn𝒫{cnwi+bn}. It can be found that the inner problem can be totally eliminated through (34) and (35), which guarantees the same value as the maximum function. Moreover, a logistic constraint that limits one efficient hyperplane in each variable region can be constructed as (36). Furthermore, constraints (37) and (38) are novelly added without loss of generality if all data points are positive. These two constraints shrink the feasible regions and accelerate the solution procedure [

29].

mincn,bni𝒟(fL(wi)-PFR(wi))2 (33)
s.t.
fL(wi)cnwi+bn    n𝒫,i𝒟 (34)
fL(wi)cnwi+bn+M(1-zin)    n𝒫,i𝒟 (35)
n𝒫zin=1    zin{0,1},iD (36)
cnwi+bn0    n𝒫,i𝒟 (37)
fL(wi)n𝒫cnwi+bn    i𝒟 (38)

where zin is a binary variable, zin=1 when the nth hyperplane takes effect for data point i; and M is a very large positive number.

As the nonlinear FRPC is a hard constraint in our planning problem, it is not necessary to force the linearized value as close to its real value as possible at all data points. In this vein, a data-driven classification idea from [

30] can be utilized to construct the linearized function. Similarly, enough fitting data points should be firstly produced with combinations of frequency-support states, whereas full enumeration is not required. The idea of this classification is applied to solve the analytical formulation, which omits the fitting errors as long as the linearized value is on the same side of the limit as its real value. On this basis, the modified DDPWL problem is proposed as (39)-(45). Specifically, constraints (40) and (41) restrict the fitting values in 𝒟1 and 𝒟3 on the same side of the limit, respectively, where the precise fitting value is not cared about. Meanwhile, all below- and above-limit data sets defined as 𝒟1, 𝒟2, and 𝒟3 can be accurately classified through (42)-(44), respectively. Note that a small number τ is introduced to guarantee the model feasibility. Based on the objective function (39), merely fitting errors are minimized in 𝒟2 [25], which not only accelerates the computation procedure, but also guarantees the conservativeness of the frequency stability.

mincn,bni𝒟2(fL(wi)-PFR(wi))2 (39)
s.t.
fL(w)<|ΔPL|    w𝒟1 (40)
fL(w)|ΔPL|    w𝒟3 (41)
𝒟1={w, PFR(w)<|ΔPL|} (42)
𝒟2={w, |ΔPL|PFR(w)<|ΔPL|+τ} (43)
𝒟3={w, |ΔPL|+τPFR(w)} (44)
(34)-(38) (45)

After deriving all PWL functions, we should deal with the maximum problem among all hyperplanes. On this basis, the big-M method [

31] can be utilized to select the largest plane in all regions subsequently. The general formulation of n hyperplanes can be expressed as (46). Finally, the constraint (29) can be transformed into a linear term |ΔPL|Ln-1. Here, u1,u2,,un-1 are binary variables to indicate which hyperplane has the maximal function value.

c1w+b1L1c1w+b1+Mu1c2w+b2L1c2w+b2+M(1-u1)L1L2L1+Mu2c3w+b3L2c3w+b3+M(1-u2)                      Ln-2Ln-1Ln-2+Mun-1cnw+bnLn-1cnw+bn+M(1-un-1) (46)

IV. Frequency-constrained Planning Model

A. Model Formulation

This section proposes a frequency-constrained GEP&ESI model, which can satisfy the required demand increase and supply enough frequency response power simultaneously. Since unit states impact system inertia and PFR directly, the short-term UC considering RPSs is inserted into the long-term planning model, which schedules sufficient frequency-support sources and promotes renewable consumption at the same time. Specifically, the objective function of the proposed model is to minimize the total cost as formulated in (47). The total cost comprises the investment costs of thermal units, wind farms, and BESS, fuel cost, unit startup cost, wind curtailment cost, load curtailment cost, load shifting-in cost, and load shifting-out cost, respectively. In this paper, the curtailable and shiftable demands are deemed as two efficient ways to obtain demand responses (DRs).

ming𝒢CCgIxg+wCwIxw+e𝒳CeIxe+s𝒮πsTht𝒯g𝒢CgGpg,t,s+g𝒢CgSUyg,t,s+wCwRpw,t,sc+b(CbDpb,t,sDC+CbFDpb,t,sFD++CbFDpb,t,sFD-) (47)

where 𝒢C, 𝒯, and 𝒮 are the sets of candidate thermal units, time intervals, and stochastic scenarios, respectively; CgI, CwI, CeI, CgG, CgSU, CwR, CbD, and CbFD are the coefficients of annual investment costs of thermal units, wind farms, and BESS, costs of fuel consumption and unit startup, scheduling of wind curtailment, load interruption, and load shifting, respectively; xg, xw, xe, pg,t,s,yg,t,s, pw,t,sc, pb,t,sDC, pb,t,sFD+, and pb,t,sFD- are decision variables including binary construction indicators of candidate thermal units, wind farms, and BESS, power output of conventional units (both thermal and hydro units), binary indicators of unit startup, curtailable power of wind farm and interrupted load, and shifting-in and shifting-out power of response demand, respectively; and πs and Th are the possibility of scenario s and the total hours of the target planning year, respectively.

This problem should subject to the following constraints.

1) Nodal Power Balance

g𝒢bpg,t,s+lto(l)=bfl,t,s-lfrom(l)=bfl,t,s+wb(pw,t,sf-pw,t,sc)+e𝒳b(re,t,sd-re,t,sc)=db,t,sDS    b,t𝒯,s𝒮 (48)

where sets 𝒢b, b, 𝒳b, and are the indices of synchronous units (both thermal and hydro units), wind farms, BESS connected to bus b, and all transmission lines, respectively; and fl,t,s and db,t,sDS are transmission power flow and the load demand after DRs, respectively.

2) DC Power Flow Constraints

DC power flow for transmission lines can be expressed as (49). Also, the maximum power flow should be limited within the line capacity by constraint (50).

fl,t,s=Bl(θo(l),t,s-θr(l),t,s)    l,t𝒯,s𝒮 (49)
Flminfl,t,sFlmax    l,t𝒯,s𝒮 (50)

where θo(l),t,s and θr(l),t,s are the voltage phase angles of origin and receive buses for line l, respectively; and Bl, Flmin, and Flmax are the susceptance, minimum and maximum power flows of line l, respectively.

3) UC Constraints

Constraint (51) implements the logistic correlation between on/off states, startup, and shutdown actions. Constraint (52) declares the exclusive correlation between startup and shutdown actions. The minimum on/off time maintained for conventional units is limited by constraints (53) and (54), respectively. Moreover, constraint (55) imposes the limit on active output power. The ramp-up and ramp-down rates for conventional generators are restricted in (56).

yg,t,s-hg,t,s=vg,t,s-vg,t-1,s    g𝒢,t𝒯,s𝒮 (51)
yg,t,s+hg,t,s1    g𝒢,t𝒯,s𝒮 (52)
k=t-UTg+1tyg,k,svg,t,s    g𝒢,t[Lg,u,NT],s𝒮 (53)
k=t-DTg+1thg,k,s1-vg,t,s    g𝒢,t[Lg,d,NT],s𝒮 (54)
Pgminvg,t,spg,t,sPgmaxvg,t,s    g𝒢,t𝒯,s𝒮 (55)
-RDgΔTpg,t,s-pg,t-1,sRUgΔT    g𝒢,t𝒯,s𝒮 (56)

where hg,t,s is the indicator of unit shutdown; vg,t,s is the indicator of on/off state of unit g; RUg and RDg are the maximum ramp-up and ramp-down power of unit g, respectively; UTg and DTg are the minimum on/off time of unit g, respectively; Lg,u and Lg,d are the duration hours of on/off states at the begging of typical days, respectively; Pgmin and Pgmax are the minimum and maximum output power offered by unit g, respectively; and NT and ΔT are the number of time intervals and the duration hour of adjacent time, respectively.

4) Wind Farm Constraints

Nowadays, in order to reduce emissions and promote investment in renewable generation, new integration incentives such as RPSs, feed-in tariffs (FITs), and production tax credits have been implemented in more than 150 countries [

32]. In this paper, for the satisfaction of the RPSs, the minimum consumption percentage γw occupied in total demand is aforehand settled in constraint (57). Also, the restriction of curtailment power with coefficient γe is also taken into account in constraint (58).

s𝒮t𝒯w(pw,t,sf-pw,t,sc)ΔTγws𝒮t𝒯bdb,t,sΔT (57)
s𝒮t𝒯wpw,t,scΔTγes𝒮t𝒯wpw,t,sfΔT (58)

where db,t,s represents the forecasted demand connected to bus b at time t in scenario s.

5) Energy Storage Constraints

The operation constraints of BESS are listed in (59)-(65). In general, the charging and discharging processes of BESS can be described as (59). Constraint (60) represents the lower and upper bounds of SOC [

33]. Besides, constraints (61) and (62) impose limitations on charging and discharging power, respectively. Constraint (63) restricts the mutual operation state. Constraint (64) stipulates that the remaining energy of BESS at the initial and end time should keep the same on representative days. Constraint (65) refers to the limitation of charging and discharging states with respect to the construction procedure. Constraint (66) proposes the preset power-headroom restrictions for BESS, which is specifically illustrated in Section II-B.

SOCe,t,s=SOCe,t-1,s+(re,t,scηec-re,t,sd/ηed)ΔT    e𝒳,t𝒯,s𝒮 (59)
SOCeminxeSOCe,t,sSOCemaxxe    e𝒳,t𝒯,s𝒮 (60)
0re,t,scRec,maxve,t,sc    e𝒳,t𝒯,s𝒮 (61)
0re,t,sdRed,maxve,t,sd    e𝒳,t𝒯,s𝒮 (62)
ve,t,sc+ve,t,sd1    e𝒳,t𝒯,s𝒮 (63)
SOCe,T0,s=SOCe,TN,s    e𝒳,t𝒯,s𝒮 (64)
ve,t,scxe,ve,t,sdxe    e𝒳,t𝒯,s𝒮 (65)
(5)-(10) (66)

where SOCe,T0,s and SOCe,TN,s are the remaining energy of BESS e at the initial and end time on representative days, respectively.

6) Cascading Hydro Generator Constraints

For hydro units, the limitations of output power, minimum on/off time, and ramp-up/down power in normal operation can be similarly constructed as thermal units. Besides, the unique constraints for cascading hydro unit can be listed as (67)-(70) [

34]. Specifically, constraint (67) restricts the water balance for hydro unit d. The water discharge and volume limits are formulated as (68) and (69), respectively. Constraint (70) restricts the initial and terminal water volume on representative days.

Vd,t,s=Vd,t-1,s+(Rd+Qd-1,t,s-Qd,t,s)    d𝒢H,t𝒯,s𝒮 (67)
vd,t,sQd,minQd,t,svd,t,sQd,max    d𝒢H,t𝒯,s𝒮 (68)
Vd,minVd,t,sVd,max    d𝒢H,t𝒯,s𝒮 (69)
Vd,0=Vd0Vd,NT=VdNT    d𝒢H (70)

where Vd,t,s and Qd,t,s are the water discharge and reservoir volume of hydro unit d, respectively; vd,t,s is the indicator of on/off state; Rd is the natural inflow to reservoir of hydro unit d; Qd,min and Qd,max are the minimum and maximum water discharges of hydro unit d, respectively;Vd,min and Vd,max are the minimum and maximum reservoir volumes of hydro unit d, respectively; and Vd0 and VdNT are the initial and terminal reservoir volumes of hydro unit d, respectively.

The water-to-power conversion of cascaded hydro units is expressed by a head-dependent function in (71) [

34]. It can be found that the output power of hydro unit d pd,t,s is solely decided on the water head level Hd,t,s, water discharge Qd,t,s, and conversion coefficient αd. Meanwhile, the water head level satisfies the correlation Hd,t,s=hd,0+βdVd,t,s, where constants hd,0 and βd are mainly determined by the physical size of reservoirs. On this basis, the output power is finally expressed as a nonlinear expression (72).

pd,t,s=αdQd,t,sHd,t,s    d𝒢H,t𝒯,s𝒮 (71)
pd,t,s=αdQd,t,s(hd,0+βdVd,t,s)    d𝒢H,t𝒯,s𝒮 (72)

7) DR Constraints

In this paper, the curtailable and shiftable demands are employed to obtain DRs [

35]. Constraint (73) restricts the curtailed power at bus b through the maximum curtailment rate λbDC. Similarly, constraint (74) models the actual demand shifting limitation through the maximum shifting rate λbFD. In the case of demand shift services, the payback effect is also considered with the payback coefficient δ, which requires the demand to be shifted out/in at bus b to satisfy the correlation in (75). In total, the relationship between the forecasted demand and DRs can be expressed as (76) and (77).

0pb,t,sDCλbDCdb,t,s    b,t𝒯,s𝒮 (73)
0pb,t,sFD+λbFDdb,t,s0pb,t,sFD-λbFDdb,t,s    b,t𝒯,s𝒮 (74)
δt𝒯pb,t,sFD-=t𝒯pb,t,sFD+    b,s𝒮 (75)
pb,t,sDR=pb,t,sDC+pb,t,sFD--pb,t,sFD+    b,s𝒮 (76)
db,t,sDS=db,t,s-pb,t,sDR    b,s𝒮 (77)

8) Frequency Constraints

Constraints (46) and (78) are utilized to guarantee enough frequency response power, while the constraint (79) requires enough system inertia to limit ROCOF. The equivalent frequency parameters can be calculated through (12), (17), and (19). Note that all equivalent frequency parameters just depend on the frequency-support states if settling a fixed unbalanced power. In this vein, the inner maximum problem can be eliminated by big-M method as introduced in Section III. Furthermore, constraints (80)-(83) preset enough power-headroom to guarantee PFR from thermal units, hydro units, wind farms, and BESS, respectively. Here, Δfss,max represents the maximum frequency deviation at the quasi-steady state, which is generally utilized to estimate the required power-headroom for PFR. If settling t in (23), Δfss,max can be estimated by (84), which mainly correlates to the values of ΔPL, SB,sys, DE, and RTmax.

ΔPL,t,sLn-1,t,s    t𝒯,s𝒮 (78)
HE,t,sfbΔPL,t,s2ROCOFmax    t𝒯,s𝒮 (79)
Pg,max-pg,t,sug,t,sKgRgPg,maxΔfss,max    g𝒢T,t𝒯,s𝒮 (80)
Pd,max-pd,t,sud,t,s1RdPd,maxΔfss,max    d𝒢H,t𝒯,s𝒮 (81)
pw,t,scuw,t,s1Rwpw,t,sfΔfss,max    w,t𝒯,s𝒮 (82)
Ee,t,sue,t,s1Re'RemaxΔfss,max    e𝒳,t𝒯,s𝒮 (83)
Δfss,max=RTmaxΔPLSB,sys(DERTmax+1) (84)

9) Binary Variable Constraints

Constraints (85)-(87) indicate that the unit states and primary frequency indicators are restricted by construction procedures. Besides, constraint (88) declares all binary variables in this model.

ug,t,svg,t,sxg    g𝒢,t𝒯,s𝒮 (85)
uw,t,sxw    w,t𝒯,s𝒮 (86)
ue,t,sxe    e𝒳,t𝒯,s𝒮 (87)
xg,xw,xe,vg,t,s,yg,t,s,hg,t,s,ve,t,sc,ve,t,sd,ug,t,s,ud,t,s,uw,t,s,ue,t,s{0,1}    g𝒢T,d𝒢H,w,e𝒳,t𝒯,s𝒮 (88)

B. Model Linearization

Unfortunately, the bilinear constraints render the proposed model unable to be dealt with by the off-the-shelf solvers. These constraints can be classified into two categories, i.e., ① the multiplier of binary and continuous variables such as (5), (6), and (9); ② the multiplier of two continuous variables such as (72). With respect to the former category, auxiliary variables denoted as qe,t,sd=re,t,sdve,t,sd and qe,t,sc=re,t,scve,t,sc are adopted to linearize constraints (5) and (6). The related auxiliary constraints are proposed as (89) and (90) through big-M method [

31]. On this basis, constraints (5) and (6) can be ultimately substituted by linear expressions (91) and (92). Similarly, constraint (9) can also be successfully linearized through the same method, which is not introduced here owing to the limitation of scope.

-M(1-ve,t,sd)qe,t,sd-re,t,sdM(1-ve,t,sd)-Mve,t,sdqe,t,sdMve,t,sd (89)
-M(1-ve,t,sc)qe,t,sc-re,t,scM(1-ve,t,sc)-Mve,t,scqe,t,scMve,t,sc (90)
0Ee,t,sdRed,maxve,t,sd-qe,t,sd    e𝒳,t𝒯,s𝒮 (91)
0Ee,t,scRed,maxve,t,sc+qe,t,sc    e𝒳,t𝒯,s𝒮 (92)

With respect to the bilinear problem with the multiplier of two continuous variables, the RLT has been used to deal with the nonlinear part in (72). This method is widely used to solve continuous factorable nonconvex optimization problems through tightly effective relaxations in a higher-dimensional space [

36]. In theory, the applied terms should satisfy the following two assumptions.

1) The nonlinear term must be expressed as the form of bilinear products.

2) All variables are bounded within the preset ranges. It can be observed from constraints (68), (69), and (72) that the nonlinear term Qd,t,sVd,t,s entirely conforms to the above conditions.

On this basis, we reformulate the nonlinear terms through substituting Qd,t,sVd,t,s by auxiliary, nonnegative, and continuous variables Ld,t,s. Meanwhile, McCormick envelopes [

36] are also used to bound these auxiliary variables as shown in (93). Finally, operation constraints of cascading hydro units can be formulated as (67)-(70), (93), and (94).

Ld,t,s-Vd,minQd,t,s-Qd,minVd,t,s+Vd,minQd,min0-Ld,t,s+Vd,maxQd,t,s-Qd,minVd,max+Qd,minVd,t,s0-Ld,t,s+Qd,maxVd,t,s-Qd,maxVd,min+Vd,minQd,t,s0Ld,t,s+Qd,maxVd,max-Qd,maxVd,t,s-Vd,maxQd,t,s0 (93)
pd,t,s=αdQd,t,shd,0+αdβdLd,t,s    d𝒢H,t𝒯,s𝒮 (94)

V. Case Study

In this section, the proposed frequency-constrained GEP&ESI model is implemented on a modified IEEE RTS-79 test system [

37]. All simulations are performed on a PC with AMD R5-3600 CPU and 16 GB RAM.

A. Test System Parameters

The modified IEEE RTS-79 test system consists of 24 buses, 33 generation units, 38 lines, and 17 load points. The demand has enlarged to be 1.1 times the original value for highlighting the needs of generation expansion. The schematic diagram of this test system is depicted in Fig. 4. All parameters of candidate thermal units can be found in Table I, where OC and IC represent operation cost and annual investment cost, respectively. Three wind farm groups named A, B, and C can be invested during the planning horizon [

38]. Specifically, group A contains three 300 MW wind farms located at buses 1, 2, and 3; group B contains three 400 MW wind farms located at buses 6, 7, and 8; group C contains three 500 MW wind farms located at buses 20, 21 and 22. All wind farms are installed with a cost of 0.5  M$/(MW·year) [39]. The parameters of candidate BESS are listed in Table II. The operation parameters of cascading hydro units can be found in [40]. Besides, the frequency parameters pertaining to thermal units, wind farms, and BESS are listed in Table III. As illustrated in Section III, the identical time constant is set as 8 s for all frequency-support devices.

Fig. 4 Schematic diagram of IEEE RTS-79 bus test system.

TABLE I Parameters of Candidate Units
Bus No.

Capacity

(MW)

NumberOC ($/MWh)

IC

($/MW)

Pg,max

(MW)

Pg,min

(MW)

1 76 2 49 200000 76 15
7 197 1 48 240000 197 38
13 197 2 48 240000 197 38
15 155 2 54 200000 155 30
16 76 2 49 190000 76 15
20 155 2 46 220000 155 20
21 375 1 48 400000 375 80
23 197 3 48 240000 197 37
TABLE II Parameters of Candidate Energy Storage
Bus No.IC ($/kWh)

Energy

(MWh)

Capacity

(MW)

Charging rateDischarging rate
1, 2, 3, 13, 16, 17, 23 500 200 100 0.9 0.875
TABLE III Frequency Parameters of Units and BESS
Units and BESSH (s)Fh (p.u.)R (p.u.)Dw (p.u.)Ke (p.u.)
Thermal unit U76 2.3 0.33 0.033
U155 4.6 0.30 0.050
U197 6.0 0.30 0.033
U375 8.0 0.35 0.050
Wind farm 3.0 15
BESS 0.050 1

Note:   U76 represents the thermal units with capacity of 76 MW.

The scenario technique is used to describe uncertainties through K-means clustering [

41]. On this basis, the historical data (8760 hourly datasets) has been clustered into four typical days to reflect both fluctuated wind power and load demand, where the forecasted data points are depicted in Fig. 5.

Fig. 5 Power proportion on four typical days.

Besides, the minimum consumption rate γw and the maximum curtailment rate γe are set as 15% and 30%, respectively, in the following case studies except otherwise stipulated. In this paper, DRs can be implemented in buses 3, 10, 13, 15, and 18, where the maximum curtailable and shiftable percentages are settled as 10%. Meanwhile, 80 $/MWh and 5 $/MWh should be paid for such curtailable and shiftable loads, respectively [

42]. The payback considers a 10% penalization in terms of energy consumption, i.e., δ=1.1. Moreover, the punishment cost of wind curtailment is settled as 150 $/MWh. In terms of frequency limits, it is considered to apply 375 MW as the instant imbalance power. The maximum frequency deviation Δfmax is limited below 0.4 Hz and ROCOFmax equals 0.5 Hz/s. The nominal frequency fb is set as 50 Hz.

B. Validation of Linearized Frequency Constraint

In this subsection, the performance of the linearized FRPC compared with the original nonlinear constraint has been examined. Two types of errors are defined as: ① Type I, satisfying the FRPC according to the linearized value but actually violating the FRPC; ② Type II, violating the FRPC according to the linearized value but actually satisfying the FRPC. The averaged fitting error errf can be calculated as:

errf=1NeωεPf(ω)-Pf,0(ω)Pf,0(ω)×100% (95)

where Ne counts the number of error points; Pf(ω) and Pf,0(ω) are the calculated values of frequency response power through linearized FRPC and initial nonlinear formulation (29), respectively; and ε is the set of error points.

Furthermore, the proposed DDPWL method is compared with three other methods that are widely used in dealing with nonlinear frequency constraints. The performances of different methods are compared in terms of the linearization error, fitting data, computation time, and the total cost. It should be noted that the bound extraction is already introduced in Section III that enumerates to seek for the worst-case frequency parameters. Comparison results of four widely-used methods can be found in Table IV. Note that PWL1 refers to the traditional PWL method in [

28] that minimizes the fitting errors at all data points. PWL2 refers to an over-conservative method, which compulsively requires the fitting value to be less than the real value. The detailed fitting model of PWL2 can be found in Appendix A. After constructing FRPC, two million random data points are produced to calculate the fitting errors of different methods.

TABLE IV Comparison Results of Four Widely-used Methods
Method

Fitting

data

Type ⅠType ⅡTotal cost (M$)

Time

(s)

NeErrf (%)NeErrf (%)
Bound extraction 2×106 0 0 0 0 390.8 9016
PWL1 4×102 4210 1.4 12100 0.90 386.4 6293
PWL2 5×104 0 0 42995 2.10 410.2 414
DDPWL 5×104 2453 0.7 381 0.02 390.6 22

In theory, no error exists for the bound extraction since the precise bound with related frequency parameters can be found by enumerating all possible data points. However, it is almost impossible to calculate the frequency response power with consideration of frequency-support states for all devices, e.g., more than 252 values exist in this test system. For the sake of simplification, two million random data points are adopted to determine the worst-case bound with more than 9016 s, which is further applied to our planning model to derive a standard total cost. As addressed in [

16], PWL1 easily causes computation obstacles since the second-order objective function minimizes the fitting errors at all data points. As simulated, only 400 data points can be successfully solved within an acceptable time (about 6293 s). Type I and Type II errors are both existed for PWL1 due to the limited information, which causes a larger possibility of insecure operation points (about 1.4%). Compared with PWL1, the computation efficiency of PWL2 is evidently improved due to the least first-order objective function. Meanwhile, PWL2 can perfectly eliminate the Type I error since all fitting values are smaller than the real values. However, an averaged Type II error of 2.1% and a total cost increment of 4.7% demonstrate its over-conservativeness under some conditions.

The above challenges can be solved by the proposed DDPWL method, which is based on the idea of data classification. It can be observed that the simulation results are promising with acceptable Type I error (0.7%) and negligible Type II error (0.02%), which almost have no impact on the total cost. Meanwhile, the computation time of DDPWL is much shorter than that of the other three methods because it only minimizes the error around the limit. In conclusion, the proposed DDPWL method can be utilized to construct the linearized FRPC with high accuracy and enhanced computation efficiency for the rest of the study.

C. Planning Results

In this sub-section, the planning results with/without frequency constraints (denoted as WFC and WOFC, respectively) are listed in Table V. Frequency responses from SGs, wind farms, and BESS are applied in WFC case. It can be found that integrating frequency constraints into the planning model does restrict the system frequency regulation. It should be emphasized that merely wind farms are expanded in WOFC case, whereas five SGs and two BESS are additionally invested in WFC case. It indicates that frequency limits facilitate the larger-scale allocation of frequency-support devices, which leads to an incremental investment cost (about 37%). Also, inserted frequency constraints bring a higher total cost (about 18.6%) to maintain frequency stability.

TABLE V Comparison of Expansion Schemes in WFC and WFOC Cases
CaseResponse typeExpansion schemesInvestment cost (M$)Fuel cost (M$)Wind curtailment cost (M$)Unit startup cost (M$)DR cost (M$)Total cost (M$)
SGsWind farmsBESS
WFC SGs+Wind farms+BESS 1(1,76), 13(1,197),21(1,375), 23(2,197) 1(300), 2(300), 3(300), 6(400), 8(400) 13, 23 102.1 281.1 5.8 1.5 0.12 390.62
WOFC None None 2(300), 6(400),8(400), 20(500) None 64.0 253.1 0.6 0 0.07 317.77

Note:   1(1,76) represents installing one thermal unit with capacity of 76 MW at bus 1; 1(300) represents installing one wind farm with capacity of 300 MW at bus 1. The scheme of BESS only lists the integration bus number.

Besides, more conventional units have to be started up and keep the minimum output power in WFC case, which increases the fuel cost from 253.1 M$ to 281.1 M$. In order to explore the effects of frequency constraints on operation schedules, the output power profiles of both existing and candidate units in all time intervals are depicted in Fig. 6. It can be observed that the online capacity of SGs in WFC case seems larger than that in WOFC case in all time intervals, while the maximum value increases from 3500 MW to 4596 MW in WFC case. The increased online capacity for SGs not only supplies inertia to limit ROCOF but also activates PFR. In addition, more power produced from SGs reduces the consumption of wind power. Also, wind farms must operate at the sub-optimal point, which causes inevitable wind curtailment to guarantee sufficient PFR. In this vein, serious wind curtailment takes place in WFC case, where an extra 5.8  M$ punishment cost is forced. The renewable curtailment strongly addresses the importance of making a trade-off between the wind consumption and frequency-support. Furthermore, the BESS will consume the spillage wind power and release energy in emergency conditions so that BESS is required to be expanded in WFC case. On the contrary, no BESS is installed in WOFC case since a smaller capacity of SGs maintains.

Fig. 6 Power profiles in all time intervals in WFC and WOFC cases. (a) WFC. (b) WOFC.

The total frequency response power from SGs, wind farms, and BESS on four representative days is shown in Fig. 7. As expected, the frequency constraints (77) and (78) would take effect to guarantee enough system frequency response power in all time intervals. On the contrary, the operation schedules in WOFC case cannot satisfy the requirement since insufficient frequency support sources are invested or kept online, which leads to evident absence power between the realistic and settled unbalanced power. As the maximum absence power arises at hour 1 on the third day, we compare the related frequency trajectories for two cases in Fig. 8. It can be found that the maximum frequency deviations for two nadir points reach 0.36 Hz and 0.86 Hz, respectively.

Fig. 7 System frequency response power in WFC and WOFC cases.

Fig. 8 Frequency trajectories of two cases in case of demand increase.

A serious frequency drop takes place in WOFC case and exceeds the permitted frequency limitation ultimately. From the perspective of ROCOF, the required inertia for this test system is 18750  MW·s for all time intervals. We can find that the minimum system inertia for WFC and WOFC cases are allocated with 22121  MW·s and 20735  MW·s, respectively, which addresses slight effects of ROCOF restriction in this test system. The above results emphasize the great significance of integrating frequency constraints into planning problems, especially for future environments with high penetration of renewable energy. It also powerfully demonstrates the effectiveness of our frequency-constrained GEP&ESI model in guaranteeing frequency stability.

D. Comparison of Different Frequency Response Modes

In this subsection, the planning results for different response modes are compared to emphasize the superiority of our proposed MSFR model. Here, four response modes are declared as: ① only SGs; ② both SGs and wind farms; ③ both SGs and BESS; ④ synthetical responses (proposed mode). The planning results are listed in Table VI. It can be observed that the proposed mode ④ derives the least total cost, where about 18  M$ (4.4%) cost savings can be obtained compared with mode ①. That addresses the great importance of taking full usage of virtual frequency responses. Meanwhile, frequency responses from wind farms seem to be more efficient than BESS when comparing modes ② and ③. It can be explained by the fact that large-scale wind farms are forced to be installed due to the consumption requirement. Thus, applying frequency response ability into wind farm control loops makes up for the reduced inertia as well as PFR from SGs, which weakens the effects of frequency limits. Also, fast PFR from BESS incurs 1.3 M$ cost savings compared with mode ①. However, normal operation restrictions affect the amounts of available response power from BESS, which evidently limits its PFR as introduced in Section II-B.

TABLE VI Planning Results with Different Response Modes
ModeIC (M$)Fuel cost (M$)Unit start-up cost (M$)Wind curtailment cost (M$)DR cost (M$)Total cost (M$)Curtailment rate (%)
SGsWind farmsBESS
35.4 64 0 305.3 3.8 0 0.04 408.54 0
26.4 64 0 295.9 1.8 3.7 0.10 391.90 0.76
32.6 64 10.4 297.2 4.8 0 0.07 409.07 0
23.7 68 10.4 281.1 1.5 5.8 0.12 390.62 1.10

E. Effects of Consumption Penetration on Planning Results

Figure 9 studies the impact of consumption penetration γw on planning results with respect to frequency constraints. Here, WWF and WOWF refer to cases with and without frequency responses from wind farms.

Fig. 9 Total cost and curtailment rate with different power penetration rates.

As can be observed from Fig. 9, the total cost of both cases gradually increases when the penetration rate ranges from 0% to 20%, since the uneconomic investment of wind farms is forced in the planning model. The cost difference between WWF and WOWF increases from 0 to 20.4 M$. The reason is that the small capacity has negligible frequency responses, which incurs no difference in a low-penetration environment. However, larger-scale wind farms enable to provide more frequency responses so that the advantage of virtual frequency response will be more evident. Besides, frequency-support incurs more curtailment power in WWF, where the curtailment rate rises initially but tends to be stable subsequently. This phenomenon addresses the restrictions of frequency-support by the installed capacity with low penetration. Since the total cost reduces to a great extent, frequency-support from virtual frequency responses can bring more cost-effective integration of wind farms when taking frequency limits into consideration.

In our planning model, the scenario technique is utilized to describe both the uncertainties of wind power and load demand by clustering the actual data set from 365 days into 4 typical days. However, the operation effectiveness that accommodates different levels of wind power and load fluctuation cannot be strictly guaranteed. In this vein, all planning schemes under different penetration rates are fixed as the input data. Then, the optimal operation simulation that considers fully actual data set is conducted. A higher possibility of wind curtailment and load shedding will be forced if the obtained schemes have violated certain operation constraints. Thus, we define violation rates of wind power VWC and load demand VLS as (96). Since all operation constraints in normal states can be totally satisfied through wind curtailment or load shedding, the normal operation problem is always feasible. However, the transient constraint of frequency nadir and ROCOF is possibly violated due to the limited frequency-support sources.

VWC=NWCNT×100%VLS=NLSNT×100% (96)

where NWC and NLS count the number of days with wind curtailment and load shedding; and NT is the number of all simulation days.

Table VII lists the simulation results with different consumption penetration rates. It can be found that the violation rate of wind power always maintains zero, while the wind consumption is verified to be larger than the requirement in realistic operation. This addresses the effectiveness of our planning model in facilitating wind power consumption, which obtains high-penetration renewable integration even merely basesd on four representative days. Regarding the violation rate of load demand VLS, certain static operation constraints are violated in view of nonzero violation rates. Nevertheless, it is so slight and negligible (less than 1%) that our planning schemes are still validated to have excellent operation effectiveness with respect to actual scenarios. Taking into account the dynamic frequency limits, the minimum frequency response power PFR and system inertia HE at all time intervals satisfy the required imbalance power (375 MW) and system inertia (18750 MW·s), respectively. This indicates that our planning schemes keep dynamic frequency stability since no frequency constraint will be violated within the actual data set.

Table VII Simulation Results with Different Consumption Penetration Rates
Penetration rate (%)VWC (%)VLS (%)Consumption rate (%)The minimum PFR (MW)Theminimum HE (MW·s)
0 0 0 0 0 0
5 0 0 0.4 0.7 0.8
10 7.5 14.6 18.8 23.1
15 375.1 375 375 375 375
20 22121 20681 20255 20255 19295

F. Effects of Unbalanced Power on Planning Results

In this subsection, we investigate the effects of unbalanced power on planning results. Figure 10 gives the total cost and curtailment rate in face of different unbalanced power. It can be observed that the total cost gradually increases with rising unbalanced power since larger-scale frequency response power is required. Similarly, PFR from wind farms plays an essential role in frequency stability due to more preset wind curtailment. However, it is limited for wind farms to participate in primary frequency regulation that can be reflected from a stable curtailment rate in the end. Besides, frequency constraints will not take effect when the unbalanced power is fairly low. For instance, a similar total cost is obtained when the unbalanced power ranges from 0 to 100 MW. Figure 11 shows the hourly system frequency response power with different unbalanced power. As the requirement ranges from 100 MW to 300 MW, the curves have significant differences, while more frequency response power will be activated. Thus, the effectiveness of our frequency-constrained GEP&ESI model in face of different unbalanced power can be totally demonstrated.

Fig. 10 Total cost and curtailment rate with different unbalanced power.

Fig. 11 Frequency response power with different unbalanced power.

VI. Conclusion

In this paper, a power-headroom constrained system frequency response model that incorporates thermal units, hydro units, wind farms, and BESS has been constructed to obtain synthetical frequency analysis. Afterwards, the multi-machine system has been transformed into a single generator by the parameter equivalence. Based on the DDPWL method, a linear coordination planning model of generation and battery energy storage has been presented to keep the frequency stability. Compared with the traditional generation planning problems, the proposed method guarantees adequate system inertia to limit ROCOF and supplies PFR to satisfy frequency nadir in future low-inertia power systems.

In our case studies, the proposed DDPWL method constructs linearized frequency constraints with high accuracy and enhanced computation efficiency. On this basis, we conduct a comparison of planning schemes in WFC and WOFC cases. Although the proposed method incurs more expansions as well as larger-scale online capacity, its effectiveness to satisfy frequency requirement is guaranteed. Moreover, the necessity of adopting virtual frequency responses for wind farms and BESSs is revealed. In such a way, more wind curtailment arises because of the sub-optimal de-loading mode, which addresses the importance of making a trade-off between wind power consumption and frequency-support. An excellent operation efficiency of the proposed model is also addressed through stochastic operation simulation. Furthermore, it can be concluded that an incremental unbalanced power incurs more installed devices since larger-scale frequency response power is required, while the effectiveness of our method to keep frequency stability is totally demonstrated.

In future work, efficient solution algorithms should be employed so that more representative days can be considered to describe uncertainty factors more precisely. Besides, more efficient control loops of frequency-support sources are expected to be applied for less conservative expansion schemes.

Appendix

Appendix A

Appendix A presents the detailed model of PWL2 methodology.

mincn,bni𝒟(fL(wi)-PFR(wi)) (A1)

s.t.

fL(wi)cnwi+bnn𝒫,i𝒟 (A2)
fL(wi)cnwi+bn+M(1-zin)n𝒫,i𝒟 (A3)
n𝒫zin=1zin{0,1},i𝒟 (A4)
fL(wi)PFR(wi) (A5)

References

1

C. Zhang, H. Cheng, L. Liu et al., “Coordination planning of wind farm, energy storage and transmission network with high-penetration renewable energy,” International Journal of Electrical Power & Energy Systems, vol. 120, pp. 1-12, Sept. 2020. [Baidu Scholar

2

S. Fang, Y. Wang, B. Gou et al., “Toward future green maritime transportation: an overview of seaport microgrids and all-electric ships,” IEEE Transactions on Vehicular Technology, vol. 69, no. 1, pp. 207-219, Oct. 2019. [Baidu Scholar

3

S. Fang, Y. Xu, S. Wen et al., “Data-driven robust coordination of generation and demand-side in photovoltaic integrated all-electric ship microgrids,” IEEE Transactions on Power Systems, vol. 35, no. 3, pp. 1783-1795, Nov. 2019. [Baidu Scholar

4

Z. Chu, U. Markovic, H. Gabriela et al., “Towards optimal system scheduling with synthetic inertia provision from wind turbines,” IEEE Transactions on Power Systems, vol. 33, no. 5, pp. 4056-4066, Apr. 2020. [Baidu Scholar

5

L. Wu and D. G. Infield, “Towards an assessment of power system frequency support from wind plant–modeling aggregate inertial response,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 2283-2291, Aug. 2013. [Baidu Scholar

6

H. Ye, W. Pei, and Z. Qi, “Analytical modeling of inertial and droop responses from a wind farm for short-term frequency regulation in power systems,” IEEE Transactions on Power Systems, vol. 31, no. 5, pp. 3414-3423, Sept. 2016. [Baidu Scholar

7

Z. Wu, W. Gao, T. Gao et al., “State-of-the-art review on frequency response of wind power plants in power systems,” Journal of Modern Power Systems and Clean Energy, vol. 6, no. 1, pp. 1-16, Jan. 2018. [Baidu Scholar

8

Y. Bao, Y. Li, B. Wang et al., “Demand response for frequency control of multi-area power system,” Journal of Modern Power Systems and Clean Energy, vol. 5, no. 1, pp. 20-29, Jan. 2017. [Baidu Scholar

9

S. Chen, T. Zhang, H. B. Gooi et al., “Penetration rate and effectiveness studies of aggregated BESS for frequency regulation,” IEEE Transactions on Smart Grid, vol. 7, no. 1, pp. 167-177, Jan. 2016. [Baidu Scholar

10

N. Nguyen, S. Almasabi, A. Bera et al., “Optimal power flow incorporating frequency security constraint,” IEEE Transactions on Industrial Applications, vol. 55, no. 6, pp. 6508-6516, Sept. 2019. [Baidu Scholar

11

Y. Y. Lee and R. Baldick, “A frequency-constrained stochastic economic dispatch model,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 2301-2312, Aug. 2013. [Baidu Scholar

12

F. Teng, V. Trovato, and G. Strbac, “Stochastic scheduling with inertia-dependent fast frequency response requirements,” IEEE Transactions on Power Systems, vol. 31, no. 2, pp. 1557-1566, Mar. 2016. [Baidu Scholar

13

F. Teng and G. Strbac, “Assessment of the role and value of frequency response support from wind plants,” IEEE Transactions on Sustainable Energy, vol. 7, no. 2, pp. 586-595, Apr. 2016. [Baidu Scholar

14

G. Zhang, E. Ela, and Q. Wang, “Market scheduling and pricing for primary and secondary frequency reserve,” IEEE Transactions on Power Systems, vol. 34, no. 4, pp. 2914-2924, Jul. 2019. [Baidu Scholar

15

P. Rabbanifar and S. Jadid, “Stochastic multi-objective security-constrained market-clearing considering static frequency of power system,” International Journal of Electrical Power & Energy Systems, vol. 54, pp. 465-480, Jan. 2014. [Baidu Scholar

16

M. Paturet, U. Markovic, S. Delikaraoglou et al., “Stochastic unit commitment in low-inertia grids,” IEEE Transactions on Power Systems, vol. 35, no. 5, pp. 3348-3458, Sept. 2020. [Baidu Scholar

17

L. Hao, J. Ji, D. Xie et al., “Scenario-based unit commitment optimization for power system with large-scale wind power participating in primary frequency regulation,” Journal of Modern Power Systems and Clean Energy, vol. 8, no. 6, pp. 1259-1267, Nov. 2020. [Baidu Scholar

18

W. Guo and J. Yang, “Modeling and dynamic response control for primary frequency regulation of hydro-turbine governing system with surge tank,” Renewable Energy, vol. 121, pp. 173-187, Jan. 2018. [Baidu Scholar

19

V. Knap, S. K. Chaudhary, D. Stroe et al., “Sizing of an energy storage system for grid inertial response and primary frequency reserve,” IEEE Transactions on Power Systems, vol. 31, no. 5, pp. 3447-3456, Sept. 2016. [Baidu Scholar

20

Y. Zhang, C. Zhao, W. Tang et al., “Profit-maximizing planning and control of battery energy storage systems for primary frequency control,” IEEE Transactions on Smart Grid, vol. 9, no. 2, pp. 712-723, Mar. 2018. [Baidu Scholar

21

M. Carrion, Y. Dvorkin, and H. Pandzic, “Primary frequency response in capacity expansion with energy storage,” IEEE Transactions on Power Systems, vol. 33, no. 2, pp. 1824-1835, Mar. 2018. [Baidu Scholar

22

S. C. Yan, Y. Zheng, and D. J. Hill, “Frequency constrained optimal siting and sizing of energy storage,” IEEE Access, vol. 7, pp. 91785-91798, Jul. 2019. [Baidu Scholar

23

C. Vrionis, V. Tsalavoutis, and A. Tolis, “A generation expansion planning model for integrating high shares of renewable energy: a meta-model assisted evolutionary algorithm approach,” Applied Energy, vol. 259, pp. 114085, Feb. 2020. [Baidu Scholar

24

V. Oree, S. Z. Hassen, and P. J. Fleming, “A multi-objective framework for long-term generation expansion planning with variable renewables,” Applied Energy, vol. 253, pp. 1135-1143, Jul. 2019. [Baidu Scholar

25

S. A. Rashidaee, T. Amraee, and M. Fotuhi-Firuzabad, “A linear model for dynamic generation expansion planning considering loss of load probability,” IEEE Transactions on Power Systems, vol. 33, no. 6, pp. 6924-6934, Nov. 2018. [Baidu Scholar

26

S. A. Rafiei, B. Mohammadi-ivatloo, S. Asadi et al., “Bi-level model for generation expansion planning with contract pricing of renewable energy in the presence of energy storage,” IET Renewable Power Generation, vol. 13, no. 9, pp. 1544-1553, Mar. 2019. [Baidu Scholar

27

Q. Shi, F. Li, and H. Cui, “Analytical method to aggregate multi-machine SFR model with applications in power system dynamic studies,” IEEE Transactions on Power Systems, vol. 33, no. 6, pp. 6355-6367, Nov. 2018. [Baidu Scholar

28

H. Ahmadi and H. Ghasemi, “Security-constrained unit commitment with linearized system frequency limit constraints,” IEEE Transactions on Power Systems, vol. 29, no. 4, pp. 1536-1545, Jan. 2014. [Baidu Scholar

29

A. Toriello and J. Vielma, “Fitting piecewise linear continuous functions,” European Journal of Operational Research, vol. 219, no. 1, pp. 86-95, Jul. 2012. [Baidu Scholar

30

Z. Chu and F. Teng, “Short circuit current constrained UC in high ibg-penetrated power systems,” IEEE Transactions on Power Systems. doi: 10.1109/TPWRS.2021.3053074. [Baidu Scholar

31

C. A Floudas, Nonlinear and Mixed-integer Optimization: Fundamentals and Applications, New York: Oxford University Press, 1995. [Baidu Scholar

32

L. Gacitua, P. Gallegos, and R. Henriquez-Auba, “A comprehensive review on expansion planning: models and tools for energy policy analysis,” Renewable Sustainable Energy Reviews, vol. 98, pp. 346-360, Sept. 2018. [Baidu Scholar

33

S. Fang, Y. Xu, Z. Li et al., “Two-step multi-objective management of hybrid energy storage system in all-electric ship microgrids,” IEEE Transactions on Vehicular Technology, vol. 68, no. 4, pp. 3361-3373, Feb. 2019. [Baidu Scholar

34

Y. Yin, T. Liu, and C. He, “Day-ahead stochastic coordinated scheduling for thermal-hydro-wind-photovoltaic systems,” Energy, vol. 187, pp. 115944, Aug. 2019. [Baidu Scholar

35

Y. Yin, T. Liu, L. Wu et al., “Frequency-constrained multi-source power system scheduling against [Baidu Scholar

contingency and renewable uncertainty,” Energy, vol. 216, pp. 119296, Feb. 2021. [Baidu Scholar

36

H. D. Sherali and W. P. Adams, A Reformulation-Linearization Technique for Solving Discrete and Continuous Nonconvex Problems, Berlin: Springer, 1998. [Baidu Scholar

37

C. Grigg, P. Albrecht, and M. Allan, “The IEEE reliability test system-1996,” IEEE Transactions on Power Systems, vol. 14, no. 3, pp. 1010-1020, Sept. 1999. [Baidu Scholar

38

D. Liu, S. Zhang, H. Cheng et al., “Accommodating uncertain wind power investment and coal-fired unit retirement by robust energy storage system planning,” CSEE Journal of Power and Energy Systems, doi: 10.17775/CSEEJPES.2019.01890 [Baidu Scholar

39

K. Tian, W. Sun, D. Han et al., “Coordinated planning with predetermined renewable energy generation targets using extended two-stage robust optimization,” IEEE Access, vol. 8, pp. 2395-2407, Dec. 2020. [Baidu Scholar

40

Y. Yin, T. Liu, and C. He, “Day-ahead stochastic coordinated scheduling for thermal-hydrowind-photovoltaic systems,” Energy, vol. 187, p. 115944, Aug. 2019. [Baidu Scholar

41

A. Soroudi and T. Amraee, “Decision making under uncertainty in energy systems: state of the art,” Renewable and Sustainable Energy Reviews, vol. 28, pp. 376-384, Aug. 2013. [Baidu Scholar

42

D. Alvarado, A. Moreira, R. Moreno et al., “Transmission network investment with distributed energy resources and distributionally robust security,” IEEE Transactions on Power Systems, vol. 34, no. 6, pp. 376-384, Nov. 2013. [Baidu Scholar