Stochastic Unit Commitment in Low-Inertia Grids

TL;DR

This paper develops a stochastic unit commitment model for low-inertia power grids, incorporating nonlinear frequency constraints and proposing an efficient linearization method, demonstrating significant impacts on operational decisions and costs.

Contribution

It introduces a novel approach to linearize nonlinear frequency nadir constraints in stochastic UC models for low-inertia grids, improving computational efficiency.

Findings

Inertia constraints significantly affect UC decisions.

Including inertia constraints increases total operational costs.

The proposed linearization method outperforms existing piece-wise linearization techniques.

Abstract

In this paper, the Unit Commitment (UC) problem in a power network with low levels of rotational inertia is studied. Frequency-related constraints, namely the limitation on Rate-of-Change-of-Frequency (RoCoF), frequency nadir and steady-state frequency error, are derived from a uniform system frequency response model and included into a stochastic UC that accounts for wind power and equipment contingency uncertainties using a scenario-tree approach. In contrast to the linear RoCoF and steady-state frequency error constraints, the nadir constraint is highly nonlinear. To preserve the mixed-integer linear formulation of the stochastic UC model, we propose a computationally efficient approach that allows to recast the nadir constraint by introducing appropriate bounds on relevant decision variables of the UC model. For medium-sized networks, this method is shown to be computationally more efficient than a piece-wise linearization method adapted from the literature. Simulation results for a modified IEEE RTS-96 system revealed that the inclusion of inertia-related constraints significantly influences the UC decisions and increases total costs, as more synchronous machines are forced to be online to provide inertial response.

Tables

Table 1. TABLE I : Computational cost of the linearization methods.

Linearization technique	Computational time [s]
PWL ($m=3$, $k=4$)	70
PWL ($m=4$, $k=4$)	7200
Bound extraction	20

Table 2. TABLE II : Parameters of the thermal plants and VSM.

Type	$H_g\left[s\right]$	$K_g[p.u.]$	$F_g[p.u.]$	$R_g[p.u.]$	$D_g[p.u.]$
Nuclear	4.5	0.98	0.25	0.04	0.6
CCGT	7.0	1.1	0.15	0.01	0.6
OCGT	5.5	0.95	0.35	0.03	0.6
VSM	6.0	1.0	-	-	0.6
Droop	-	1.0	-	0.05	-

Table 3. TABLE III : Comparison of the total number of dispatched generators through UC for each hour.

Hour	65	66	67	68	69	70	71	72	73
w/o FC	6	5	4	4	4	4	4	4	4
w/ FC	6	5	10	10	10	10	10	10	10

Table 4. TABLE IV : Unit commitment costs [ $ ] delimited-[] currency-dollar \mathrm{[}\$] breakdown on day 3.

Scenario	Total costs	Start-up	Operation	Reserves
w/o FC	$\mathrm{410\hspace{0.17em}545}$	$788$	$\mathrm{405\hspace{0.17em}561}$	$\mathrm{4\hspace{0.17em}196}$
w/ FC	$\mathrm{432\hspace{0.17em}383}$	$\mathrm{2\hspace{0.17em}248}$	$\mathrm{420\hspace{0.17em}886}$	$\mathrm{9\hspace{0.17em}249}$
Difference	$5.32\%$	$185\%$	$3.78\%$	$120\%$

Equations

$$G(s)$$

$$\displaystyle+\underbrace{\sum\limits_{d\in\mathcal{J}_{d}}\dfrac{K_{d}}{R_{d}(1+sT_{d})}}_{\text{droop converters}}+\underbrace{\sum\limits_{v\in\mathcal{J}_{v}}\dfrac{sM_{v}+D_{v}}{1+sT_{v}}}_{\text{VSM converters}}\bigg{)}^{-1}.$$

$$G(s)=\frac{1}{MT}\frac{1+sT}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2},$$

$${\omega }_n=\sqrt{\frac{D+R_g}{MT}},\zeta =\frac{M+T(D+F_g)}{2\sqrt{MT(D+R_g)}},$$

$$\Delta f_{\max }$$

$${\dot{f}}_{\max }$$

$$\Delta f_{\mathrm{ss}}$$

$$\left|\frac{f_b\Delta P}{D+R_g}\left(1+\sqrt{\frac{T(R_g-F_g)}{M}}{e}^{-\zeta {\omega }_n{t}_m}\right)\right|$$

$$\left|\frac{f_b\Delta P}{M}\right|\le {\dot{f}}_{\lim },$$

$$\left|\frac{f_b\Delta P}{D+R_g}\right|\le \Delta f_{\mathrm{ss},\lim },$$

$$\displaystyle\min_{\Psi}\sum_{\eta}\bigg{(}$$

$$\displaystyle-\Delta f_{\mathrm{max}}\left(R_{g}^{(\eta)},F_{g}^{(\eta)},M^{(\eta)}\right)\bigg{)}^{2}$$

$$F_g\ge F_g^{\lim },R_g\ge R_g^{\lim },M\ge M^{\lim }.$$

$$\pi [A(1,\tau )]=\exp (-{\lambda }_1\tau )(\exp ({\lambda }_1)-1),$$

$$\pi [A(1,\tau )]=\pi [A(2,\tau )]=\dots =\pi [A(\mathcal{K},\tau )],$$

$$\pi [B(1)]=\exp (-{\lambda }_1\tau ),$$

$$\pi [B(1)]=\pi [B(2)]=\dots =\pi [B(\mathcal{K})].$$

$${\pi }_{c_0}$$

$${\pi }_{c_{\kappa }}=\pi [A(1,\tau )]\prod _{\begin{array}{c}y=1\\ y\ne k\end{array}}^{\mathcal{K}}\pi [B(y)],\forall \kappa =1,...,\mathcal{K}.$$

$$\displaystyle\min_{\Phi}\,\,\sum_{t\in\mathcal{T}}\sum_{i\in\mathcal{I}}\big{(}C^{\mathrm{SU}}_{i}y_{it}+C^{\mathrm{SD}}_{i}z_{it}+C_{i}p_{it}\big{)}\ +$$

$$\displaystyle\sum_{t\in\mathcal{T}}\sum_{\xi\in\mathcal{E}}\pi_{\xi}\Big{[}\sum_{i\in\mathcal{I}}\big{(}C^{+}_{i}r^{+}_{i\xi t}-C^{-}_{i}r^{-}_{i\xi t}\big{)}+\sum_{n\in\mathcal{N}}C^{\mathrm{sh}}\,l^{\mathrm{shed}}_{n\xi t}\Big{]}$$

$$\sum _{i\in {\mathcal{I}}_n}{p}_{it}+\sum _{j\in {\mathcal{J}}_n}{w}_{jt}-D_{nt}-$$

$$\sum _{m:(n,m)\in \mathcal{L}}{B}_{nm}({\hat{\delta }}_{nt}-{\hat{\delta }}_{mt})=0,\forall n,t,$$

$$B_{nm}({\hat{\delta }}_{nt}-{\hat{\delta }}_{mt})\le {\overline{f}}_{nm},\forall (n,m)\in \mathcal{L},t,$$

$$u_{it}-u_{i(t-1)}\le u_{i{\tau }_i^1},\forall i,t,$$

$$u_{i(t-1)}-u_{it}\le 1-u_{i{\tau }_i^0},\forall i,t,$$

$$y_{it}\ge u_{it}-u_{i(t-1)},\forall i,t,$$

$$z_{it}\ge u_{i(t-1)}-u_{it},\forall i,t,$$

$$\displaystyle\sum_{i\in\mathcal{I}_{n}}\Big{[}r^{+}_{i\xi t}-r^{-}_{i\xi t}-p_{it}(1-\alpha_{i\xi t})\Big{]}$$

$$+\sum _{m:(n,m)\in \mathcal{L}}{B}_{nm}({\hat{\delta }}_{nt}-{\tilde{\delta }}_{n\xi t}-{\hat{\delta }}_{mt}+{\tilde{\delta }}_{m\xi t})$$

$$+\sum _{j\in {\mathcal{J}}_n}(W_{j\xi t}^{\ast }-w_{jt}-w_{j\xi t}^{\mathrm{spill}})+l_{n\xi t}^{\mathrm{shed}}=0,\forall n,\xi ,t,$$

$$p_{it}+r_{i\xi t}^{+}\le \overline{P_i}{u}_{it},\forall i,\xi ,t,$$

$$p_{it}-r_{i\xi t}^{-}\ge \underline{P_i}{u}_{it},\forall i,\xi ,t,$$

$$p_{it}-p_{i(t-1)}+r_{i\xi t}^{+}-r_{i\xi (t-1)}^{+}\le R_i^{\mathrm{U}},\forall i,\xi ,t,$$

$$p_{it}-p_{i(t-1)}-r_{i\xi t}^{-}+r_{i\xi (t-1)}^{-}\ge -R_i^{\mathrm{D}},\forall i,\xi ,t,$$

$$r_{i\xi t}^{+}\le R_i^{+}{\alpha }_{i\xi t},\forall i,\xi ,t,$$

$$r_{i\xi t}^{-}\le R_i^{-}{\alpha }_{i\xi t},\forall i,\xi ,t,$$

$$B_{nm}({\tilde{\delta }}_{nt}-{\tilde{\delta }}_{mt})\le {\overline{f}}_{nm},\forall (n,m)\in \Lambda ,t,$$

$$w_{j\xi t}^{\mathrm{spill}}\le w_{j\xi t},\forall j,\xi ,t,$$

$$l_{n\xi t}^{\mathrm{shed}}\le D_{nt},\forall n,\xi ,t,$$

$$k_{i\xi t}=\frac{P_i{K}_i}{{\sum }_{i\in \mathcal{I}}{P}_i+W_v+W_d}{u}_{it}{\alpha }_{i\xi t},\forall i,\xi ,t,$$

$$F_{\xi t}=\sum _{i\in \mathcal{I}}\frac{F_i{k}_{i\xi t}}{R_i},\forall t,\xi ,$$

$$R_{\xi t}=\sum _{i\in \mathcal{I}}\frac{k_{i\xi t}}{R_i},\forall t,\xi ,$$

$$M_{\xi t}=\sum _{i\in \mathcal{I}}2H_g{k}_{i\xi t},\forall t,\xi ,$$

$$\frac{{\dot{f}}_{\lim }}{f_b}(M_{\xi t}+M_v)\ge \Delta P_{\xi t},\forall t,\xi ,$$

$$F_{\xi t}\ge F_{\xi t}^{\lim },R_{\xi t}\ge R_{\xi t}^{\lim },M_{\xi t}+M_v\ge M_{\xi t}^{\lim },\forall t,\xi ,$$

$$\frac{\Delta f_{\mathrm{ss},\lim }}{f_b}(D+R_{\xi t})\ge \Delta P_{\xi t},\forall t,\xi ,$$

$$p_{it}\ge 0,\forall i,t;w_{jt}\ge 0,\forall j,t;{\hat{\delta }}_{nt}\ge 0,\forall n,t;k_{i\xi t}\ge 0,\forall i,\xi ,t;$$

$$r_{i\xi t}^{+},r_{i\xi t}^{-}\ge 0,\forall i,\xi ,t;l_{n\xi t}^{\mathrm{shed}}\ge 0,\forall n,\xi ,t;w_{j\xi t}^{\mathrm{spill}}\ge 0,\forall j,\xi ,t;$$

Full text

Nomenclature

The main notation used in this paper is introduced below. Additional symbols are defined in the paper where needed. All symbols are augmented by index $t$ when referring to different time periods.

-A Sets and Indices

• $\ell\in\mathcal{L}$

Set of transmission lines.

• $\xi\in\mathcal{E}$

Set of scenarios $\xi=\{c,\omega\}$ including generation outages ($c$) and wind power uncertainty ($\omega$).

• $i\in\mathcal{I}$

Set of conventional generation units.

• $j\in\mathcal{J}$

Set of converter-based (i.e., wind) generation units.

• $n\in\mathcal{N}$

Set of nodes.

• $\mathcal{I}_{n}$

Set of conventional generation units located at bus $n$.

• $\mathcal{J}_{n}$

Set of converter-based units located at bus $n$.

-B Decision variables

• $\hat{\delta}_{n}$

Day-ahead voltage angle at node $n$ $[\mathrm{rad}]$.

• $\tilde{\delta}_{n\xi}$

Real-time voltage angle at node $n$ in scenario $\xi$ $[\mathrm{rad}]$.

• $F_{\xi t}$

Global fraction of total power generated by high-pressure turbines in scenario $\xi$ $[\mathrm{p.u.}]$.

• $k_{i\xi t}$

Scaled power gain factor of conventional unit $i$ $\xi$ $[\mathrm{MW}]$.

• $l_{n\xi t}^{\mathrm{shed}}$

Shedding of load at node $n$ in scenario $\xi$ $[\mathrm{MW}]$.

• $M_{\xi t}$

Global system inertia in scenario $\xi$ $[\mathrm{p.u.}]$.

• $p_{it}$

Day-ahead dispatch of conventional unit $i$ $[\mathrm{MW}]$.

• $r_{i\xi t}^{+/-}$

Up-/Downward reserve deployment of unit $i$ in scenario $\xi$ $[\mathrm{MW}]$.

• $R_{\xi t}$

Global system droop factor in scenario $\xi$ $[\mathrm{p.u.}]$.

• $u_{it}$

Commitment variable of conventional unit $i$.

• $w_{jt}$

Day-ahead dispatch of wind power unit $j$ $[\mathrm{MW}]$.

• $w_{j\xi t}^{\mathrm{spill}}$

Wind spillage of unit $j$ in scenario $\xi$ $[\mathrm{MW}]$.

• $y_{it}$

Start-up variable of conventional unit $i$ .

• $z_{it}$

Shut-down variable of conventional unit $i$ .

-C Parameters

• $\alpha_{i\xi t}$

Outage parameter of conventional unit $i$ in scenario $\xi$.

• $\Delta P_{\xi t}$

Size of power outage in scenario $\xi$ $[\mathrm{p.u.}]$.

• $\pi_{\xi}$

Probability of occurrence of scenario $\xi$.

• $B_{nm}$

Susceptance of transmission line $(n,m)$ $[\mathrm{S}]$.

• $C_{i}$

Day-ahead price offer of unit $i$ [\mathrm{\/MWh}]$.

• $C_{i}^{\mathrm{SU/SD}}$

Start-up/Shut-down price offer of unit $i$ [\mathrm{\}]$.

• $C_{i}^{\mathrm{+/-}}$

Up-/Downward reserve price offer of unit $i$ [\mathrm{\/MWh}]$.

• $C^{\mathrm{sh}}$

Value of lost load [\mathrm{\/MWh}]$.

• $D_{nt}$

Demand at node $n$ $[\mathrm{MW}]$.

• $\overline{f}_{nm}$

Capacity of transmission line $(n,m)$ $[\mathrm{MW}]$.

• $\overline{P_{i}}/\underline{P_{i}}$

Active power limits of conventional unit $i$ $[\mathrm{MW}]$.

• $R_{i}^{\mathrm{+/-}}$

Up-/Downward reserve capacity of unit $i$ $[\mathrm{MW}]$.

• $R^{\mathrm{U}/\mathrm{D}}_{i}$

Ramp up/down limits of conventional unit $i$ $[\mathrm{MW/h}]$.

• $W_{j\xi t}^{*}$

Wind power realization of unit $j$ in scenario $\xi$ $[\mathrm{MW}]$.

• $W_{d}$

Total capacity of droop control units $[\mathrm{MW}]$.

• $W_{v}$

Total capacity of VSM units $[\mathrm{MW}]$.

I Introduction

With increasing penetration of renewable energy sources, system operators face new challenges in order to ensure power grid stability. One of these challenges is frequency stability due to a loss of generation or a large variation of load. In traditional power systems, synchronous generators (e.g., hydro or steam turbines) provide rotational inertia through stored kinetic energy in their rotating mass (turbine system and rotor). This energy is important to stabilize the system as it ensures slower frequency dynamics and reduces the Rate of Change of Frequency (RoCoF) in case of a generation-demand imbalance [1]. In the future, with more generation coming from wind and solar power, the ability of the system to maintain the frequency within the acceptable range is diminished. Indeed, photovoltaic systems are connected to the grid through inverters, which do not exhibit rotational inertia. Even in the case of inverter-interfaced wind generators, the inverter electrically decouples the rotor’s rotational inertia from the system [2].

Transmission System Operators (TSOs) around the world are concerned with the stability issues associated with large penetration of renewable energy in their systems. In the United States, the Electric Reliability Council of Texas (ERCOT) has studied the effect of low inertia on the security and reliability of the grid [3]. The Irish TSO, EirGrid, is designing ancillary services to remunerate providers of rotational or synthetic inertia [4]. Furthermore, EirGrid currently imposes limits on the maximum instantaneous penetration of variable RES with respect to the total load demand at any point in time.

In systems with low rotational inertia, TSOs must impose minimum inertia requirements in order to secure frequency stability and avoid system collapse in case of a severe fault or a sudden mismatch between generation and demand. With such new requirements, the traditional Unit Commitment (UC) problem, i.e., the day-ahead scheduling process to decide which generators will be committed, may be affected as more Synchronous Generators (SGs) could be dispatched for the sole purpose of providing inertia.

Several papers have approached the problem of including inertia requirements in the UC problem. In [5, 6, 7], the authors use the swing equation of Center-of-Inertia (CoI), which allows them to derive the RoCoF constraint and study its effect on the UC schedule. However, this approach oversimplifies the problem as it neglects metrics related to frequency deviation from the setpoint. This problem was addressed in [8] and [9] with the inclusion of a constraint limiting the post-disturbance maximum frequency deviation (i.e., frequency nadir). In [8], the analytic form of frequency nadir as a function of active power disturbance is derived using a system frequency model obtained from [10]. The nadir expression is then linearized and added to the UC model, while considering a fixed sudden load increase. On the other hand, the authors of [9] bypass the explicit modeling of turbine and governor control, as well as their impact on frequency dynamics, by imposing strict assumptions on system damping and total frequency response provision at each node. Moreover, they look at the impact of wind uncertainty on inertia requirements. While the simplifications proposed in [9] enable the formulation of a nadir constraint without the explicit consideration of second-order frequency dynamics, they oversimplify the actual control implementation and disregard the aggregate impact of governor damping. Furthermore, none of the aforementioned studies incorporate the converter interface of RES and the impact of respective control schemes on the UC formulation.

This study builds on the work of [8] and [9], and improves on it in several ways. First, we improve the frequency dynamics model in [8] by including the state-of-the-art converter control schemes of inverter-based generation, more specifically Virtual Synchronous Machine (VSM) and droop control. In contrast to the existing literature, where SG inertia and damping constants are usually numerically modified in order to compensate for high RES integration, we analyze a realistic model of a low-inertia system comprising both SG and converter dynamic models. This allows us to derive detailed analytic expressions of relevant frequency metrics as functions of multiple system variables (e.g., inertia, damping, aggregate droop gain, etc.) to be determined by the UC model, as opposed to the approach in [9] where inertia constant was the only decision variable of interest. Moreover, in addition to frequency nadir and RoCoF, we incorporate the limitation on quasi steady-state frequency deviation into the UC formulation. Secondly, a more straightforward method is proposed to extract bounds for decision variables contributing to frequency nadir, which allows us to incorporate the non-linear nadir constraint in the UC problem in a more efficient way compared to [8]. Furthermore, similar to [9], this paper includes both the wind uncertainty and potential loss of generation in the UC model. However, we present a more comprehensive approach towards event probability computation and structuring of the scenario tree for the two-stage stochastic UC problem.

The rest of the paper is structured as follows. In Section II, the derivation of post-contingency frequency dynamics in a low-inertia, multi-machine system is discussed. The obtained time-domain, analytic expressions are then linearized and incorporated into a stochastic UC formulation in Section III. Subsequently, the modeling of uncertainties, namely equipment failure and wind power, in the form of probabilistic scenarios is presented in Section IV. Section V provides the mathematical formulation of the stochastic UC problem. Finally, Section VI presents and discusses the simulation results using a modified version of the IEEE RTS-96 system, whereas Section VII draws the main conclusions and discusses the outlook of the study.

II Low-Inertia System Frequency Dynamics

II-A Inertial Response and Primary Frequency Control Model

We first focus on deriving a simplified, but sufficiently accurate, uniform frequency response model of a low-inertia system previously introduced in [11]. Let us consider a system comprised of traditional $(i\in\mathcal{I})$ and converter-based $(j\in\mathcal{J})$ generators depicted in Fig. 1.

The generator dynamics are described by the swing equation, with $M_{g}$ and $D_{g}$ denoting the normalized inertia and damping constants corresponding to the synchronous generators’ CoI. The low-order model proposed in [12] is used for modeling the governor droop and turbine dynamics; $T_{i}$ are the turbine time constants, $R_{i}$ and $K_{i}$ are the respective droop and mechanical power gain factors, while $F_{i}$ refers to the fraction of total power generated by the turbines of synchronous machines. Furthermore, we incorporate the impact of grid-forming converters, as they are the only type of power electronic-interfaced units providing frequency support [13, 14]. A particular focus is set on droop $(d\in\mathcal{J}_{d}\subseteq\mathcal{J})$ and VSM $(v\in\mathcal{J}_{v}\subseteq\mathcal{J})$ control schemes, as two of the currently most prevalent emulation techniques in the literature, which in fact have equivalent properties in the grid-forming mode of operation [15]. Here, $T_{d}=T_{v}\equiv T_{j}$ are the time constants of all converters, $R_{d}$ and $K_{d}$ are the respective droop and electrical power gain factors, whereas $M_{v}$ and $D_{v}$ denote the normalized virtual inertia and damping constants of VSM converters.

II-B Analytic Derivation of Frequency Metrics

From Fig. 1 we can now derive a transfer function $G(s)$ of a general-order system dynamics, as follows:

[TABLE]

Assuming similar time constants $(T_{i}\approx T)$ of all synchronous machines, usually 2-3 orders of magnitude higher than the ones of converters, justifies the approximation $T\gg T_{j}\approx 0$. Now we can transform (1) into the following expression:

[TABLE]

where the natural frequency $(\omega_{n})$ and damping ratio $(\zeta)$ are

[TABLE]

and parameters $(M,D)$ and $(F_{g},R_{g})$ represent weighted system and synchronous generator averages, respectively. More details on mathematical formulation can be found in [11].

Assuming a stepwise disturbance in the electrical power $\Delta P_{e}(s)=-\Delta P/s$, we can derive the time-domain expression for frequency deviation as well as the nadir ($\dot{f}_{\mathrm{max}}$), RoCoF ($\dot{f}_{\mathrm{max}}$) and steady-state deviation ($\Delta f_{\mathrm{ss}}$) frequency metrics:

[TABLE]

with the introduction of the new variable $\omega_{d}=\omega_{n}\sqrt{1-\zeta^{2}}$.

The accuracy of the proposed model has already been investigated and verified in [11]. We can conclude that the frequency metrics of interest are directly dependent on the average system parameters $M$, $D$, $R_{g}$ and $F_{g}$, and thus they could be regulated through the UC model. In particular, RoCoF and steady-state deviation can be explicitly controlled via $\dot{f}_{\mathrm{max}}\sim M^{-1}$ and $\Delta f_{\mathrm{ss}}\sim(D+R_{g})^{-1}$, while nadir can be modeled using a highly non-linear function $\Delta f_{\mathrm{max}}\left(M,D,R_{g},F_{g}\right)$.

III Formulation of Frequency Constraints

The aforementioned frequency expressions in (4) are incorporated as constraints into the stochastic UC problem, converted into SI and bounded by prescribed ENTSO-e thresholds [16], as follows:

[TABLE]

with $f_{b}=50\,\mathrm{Hz}$ being the base frequency; $\Delta f_{\mathrm{lim}}=0.4\,\mathrm{Hz}$ is the Under-Frequency Load Shedding (UFLS) trigger, while $\dot{f}_{\mathrm{lim}}=0.5\,\mathrm{Hz}/\mathrm{s}$ and $\Delta f_{\mathrm{ss,lim}}=0.2\,\mathrm{Hz}$ are the maximum permissible RoCoF and steady-state frequency deviation.

Constraints (5b) and (5c) are linear, unlike the non-linear frequency nadir constraint (5a). In order to avoid the high computational burden of a Mixed-Integer Non-Linear Program (MINLP) formulation and have a measurable optimality gap, we use a linear approximation of (5a) which allows us to maintain a Mixed-Ineger Linear Program (MILP) formulation of the stochastic UC problem.

III-A Piece-wise Linearization of Nadir Expression

The study in [8] proposes a Piece-Wise Linearization (PWL) technique for obtaining a linearized expression for frequency nadir in order to subsequently integrate it into a UC problem. To improve clarity, this technique is outlined here and in Appendix A before comparing its computational burden against our proposed approach introduced in Section III-B. Let us recall from Section II-B that the frequency nadir expression is a function of four variables $(R_{g},F_{g},M,D)$, and as such too complicated to be directly handled by the PWL. Considering that the aggregate damping constant is of the form $D(D_{i},D_{v},R_{d})$, with respective damping and droop gains usually strictly prescribed within narrow ranges by the system operator, it is justifiable to assume a constant $D$ and therefore becomes $\Delta f_{\mathrm{max}}(M,R_{g},F_{g})$. Hence, the PWL formulation aims to minimize the following objective function

[TABLE]

with $\Psi=\{a_{\nu},b_{\nu},c_{\nu},d_{\nu},\forall\nu\}$, being the set of optimization variables, $\eta$ denoting the evaluation point and $\nu$ referring to the number of PWL segments. The objective function (6) penalizes the difference between the appropriate PWL segment and the nadir function at all evaluation points. Given the convex nature of the nadir function, the inner $\max$ operator chooses the appropriate PWL segment for each evaluation point by looking at which segment is closest to the curve at that specific point. To improve understanding, we provide an illustration of the PWL method from [8] on a simple one variable function in Appendix A.

Upon obtaining the optimal solution of the model in (6), denoted as $(a^{*}_{\nu},b^{*}_{\nu},c^{*}_{\nu},d^{*}_{\nu})$, the nadir constraint can be integrated into the MILP UC model by adding a set of inequalities described in Appendix A, along with the nadir threshold constraint of the form $f_{b}\,t_{3}\leq\Delta f_{\mathrm{lim}}$. The results for the approximation of frequency nadir function for a test system of 20 generators described in Section VI are shown in Fig. 2, where a loss of the largest unit is considered. Note that Fig. 2 showcases the surface plot for a fixed inertia constant $M$ and thus ignores one degree of freedom. The original surface is presented in blue, whereas its PWL-approximation segments are the planes depicted in various colors. It is important to note that the optimization problem (6) is computationally intensive and thus in order to obtain results within reasonable computational time, the number of PWL segments used for the approximation as well as the number of evaluation points have to be kept low.

III-B Extracting Bounds on Relevant Variables

An alternative approach for linearizing the nadir constraint and integrating it into the UC problem is to confine the values of $R_{g}$, $F_{g}$, $M$ and $D$ within a plausible range to guarantee that the nadir threshold in (5a) is not violated. With this approach the damping variable $D$ can easily be included and does not need to be set constant. The scatter plot presented in Fig. 3 reflects all possible values of frequency nadir after the loss of the largest generator, for the same system as in Section III-A. In the general case, for a system that comprises $|\mathcal{I}|$ generators, there will be $2^{|\mathcal{I}|-1}$ possible generator commitment combinations after a generator outage. By obtaining the set of these dispatch combinations, the values of $R_{g}$, $F_{g}$ and $M$ at which the UFLS threshold is not violated can be extracted, corresponding to the points below the shaded plane in Fig. 3. Subsequently, these values are used to substitute the nadir constraint in the unit commitment as follows:

[TABLE]

Table I provides a comparison of the proposed method to the PWL technique, in terms of computational time that is needed to obtain the equivalent linear nadir equations for a single value of $\Delta P$. It is clear that the PWL is more time intensive, especially when aiming for an increased precision. It should be noted though that for very large systems, the calculations of $2^{|\mathcal{I}|-1}$ combinations for the bound extraction method would become more computationally expensive. Indeed, the computation time increases by a factor of $2^{\Delta|\mathcal{I}|}$ for every additional $\Delta|\mathcal{I}|$ generators included in the system. For the purposes of this paper, the proposed bound extraction method will be used as it is significantly faster and introduces less error when applied to the 20-generator test system under investigation.

IV Modeling equipment-failure and wind power uncertainties

This section describes the modeling of uncertainty pertaining to equipment failure and wind power production during power system operation. The uncertain nature of wind power production is modeled using a set of scenarios $\Omega$ that captures the spatio-temporal interdependence of forecast errors, for every wind farm location and during the whole scheduling horizon. Each wind power realization scenario $\omega$ has the same probability of occurrence denoted as $\pi_{\omega}$.

In terms of equipment failure uncertainty, we consider as the set of credible contingencies, the unforeseen outages of synchronous generators, whereas transmission assets are assumed to be $100\,\%$ reliable. In order to reduce the computational burden, we follow the assumption from [17] considering that the generation outages happen at a discrete time period, while failed assets remain unavailable for the rest of the scheduling horizon, i.e., the Mean Time To Repair (MTTR) is greater than the scheduling horizon of the day-ahead electricity market. For the purpose of assessing the impact of frequency constraints on the unit commitment schedule, we consider as contingency period the one in which the power system faces the highest wind power penetration as scarcity of inertia is most likely to occur in this time due to the displacement of synchronous generators from the day-ahead schedule.

To calculate the probability $\pi_{c}$ associated with contingency scenario $c$, we index the set of credible contingencies by $\kappa=1,2,\dots,\mathcal{K}$. We denote by $A(\kappa,\tau)$ the random event of contingency $\kappa$ happening within time period $\tau$. Random event $B(\kappa)$ corresponds to contingency $\kappa$ not occurring during the entire scheduling horizon. We further denote as $\lambda_{\kappa}$ the inverse of the Mean Time To Failure $(\mathrm{MTTF})$ $\kappa$, i.e., $\lambda_{\kappa}=1/\mathrm{MTTF}$. Considering we are looking at only one hour in which the outages may occur, the probability $\pi_{c}$ for each contingency scenario $\kappa$ is derived from the probabilities of occurrence of random events $A(\kappa,\tau)$ and $B(\kappa)$ that are calculated using the following expressions according to [18]:

[TABLE]

Assuming statistical independence between all contingencies, the probability $\pi_{c_{0}}$ of the no-contingency scenario is equal to

[TABLE]

while the probability $\pi_{c_{\kappa}}$ of losing a generator is equal to

[TABLE]

It should be noted that the sum of probabilities $\pi_{c_{0}}$ and $\pi_{c_{\kappa}}$ is lower than 1, since sequential contingencies are not considered. For instance, setting $\mathrm{MTTF}$ equal to $1000\,\mathrm{h}$ for all generators, we obtain from (8) $\pi[A(\kappa,\tau)]=0.9995\!\times\!10^{-3}$ and $\pi[B(\kappa)]=0.9990$. According to (9) and (10) we obtain respectively $\pi_{c_{0}}=0.9960$ and $\pi_{c_{\kappa}}=0.9965\times 10^{-3}$ and thus $\sum_{\kappa=0}^{4}\pi_{c_{\kappa}}=0.9999\approx 1$.

Combining the scenarios modeling the equipment failure and wind power uncertainty into a single scenario set $\mathcal{E}$, we define each scenario $\xi$ as a pair of contingency $c$ and wind power realization $\omega$. For each $\xi=\{c,\omega\}$ the corresponding probability of occurrence is given as $\pi_{\xi}=\pi_{\omega}\cdot\pi_{c}$ and $\sum_{\xi\in\mathcal{E}}\pi_{\xi}\approx 1$, assuming that equipment outages and wind power production are statistically independent events. The structure of the scenario set $\mathcal{E}$ used in the stochastic UC formulation is illustrated as the scenario tree shown in Fig. 4, for $|\mathcal{K}|$ contingencies and $|\mathcal{W}|$ wind power scenarios.

V Stochastic Unit Commitment

This section provides the mathematical formulation of the stochastic unit commitment [19], with an addition of frequency-related constraints. The proposed model is a two-stage stochastic optimization problem which can be written as:

[TABLE]

where $\Phi=\{p_{it},\,u_{it},\,y_{it},\,z_{it},\,\forall i,t;\;w_{jt},\forall j,t;\;\hat{\delta}_{nt},\forall n,t;$ $\tilde{\delta}_{n\xi t},\forall n,\xi,t;\,r^{+}_{i\xi t},\,r^{-}_{i\xi t},\forall i,\xi,t;\,w^{\mathrm{spill}}_{j\xi t},\forall j,\xi,t;\,l^{\mathrm{shed}}_{n\xi t},\forall n,\xi,t;\,$ $k_{i\xi t},\forall i,\xi,t;\,F_{\xi t},\,R_{\xi t},\,M_{\xi t},\forall\xi,t\}$ is the set of optimization variables.

The objective function (11a) to be minimized is the total expected system cost that comprises the day-ahead energy and the real-time balancing costs. The day-ahead component consists of the fuel costs $C_{i}$ as well as the start-up and shut-down costs. The real-time component includes the re-dispatch cost from the deployment of upward and downward reserves based on the corresponding offer prices $C^{+}_{i}$ and $C^{-}_{i}$, as well as the involuntary load shedding at the value of lost load $C^{\mathrm{sh}}$.

Equation (11b) enforces the nodal power balance of the day-ahead schedule, while network power flows at the day-ahead stage are restricted by the transmission capacity limits in (11c). Constraints (11d)-(11e) model the minimum online and offline time of conventional units based on commitment variable $u_{it}$, where parameters $\tau^{1}_{i}$ and $\tau^{0}_{i}$ are defined as $\tau_{i}^{1}=\min\{t+T^{1}_{i}-1,\ T\}$ and $\tau_{i}^{0}=\min\{t+T^{0}_{i}-1,\ T\}$, and $T^{1}_{i}$ and $T^{0}_{i}$ denote the duration that unit $i$ should remain online and offline, respectively. Constraints (11f)-(11g) model the start-up and shut-down of conventional units using the binary variables $y_{it}$ and $z_{it}$, respectively. The real-time power balance for every uncertainty realization $\xi$ is enforced by constraint (11h). Parameter $\alpha_{i\xi t}$ models the availability of the generators to provide reserves, i.e., $\alpha_{i\xi t}$ is equal to 1 if generator $i$ at scenario $\xi$ and time $t$ is online and able to provide reserves and zero otherwise. The scheduled energy production and the deployment of upward ($r^{+}_{i\xi t}$) and downward ($r^{-}_{i\xi t}$) reserves in each scenario $\xi$ are bounded by the generation capacity limits of each unit by constraints (11i)-(11j), whereas constraints (11k)-(11l) enforce the upward and downward ramping limits accounting for the real-time reserve activation. Constraints (11m)-(11n) account for the limits of reserve capacity offers. Transmission capacity limits during real-time operation are enforced by constraint (11o), whereas wind spillage $w^{\mathrm{spill}}_{j\xi t}$ and load shedding $l^{\mathrm{shed}}_{n\xi t}$ are bounded by the wind power realization and the nodal demand through constraints (11p) and (11q), respectively.

The set of constraints (11r)-(11w) models the frequency limits of the power system. The equality constraint (11r) defines $k_{i\xi t}$ as the gain factor $K_{i}$ of generator $i$ scaled by the ratio of its capacity over the total system capacity, which in turn is multiplied by the binary variable $u_{it}$ and the parameter $\alpha_{i\xi t}$ to indicate that a unit can only provide inertial response if it is committed and does not face an outage. Similarly, constraints (11s)-(11u) define average system variables for power fraction, droop and inertia, respectively. Constraint (11v) enforces the RoCoF limit, while nadir equivalent and quasi steady-state frequency bounds are imposed by constraints (11w) and (11x). Finally, constraints (11y) are variable declarations.

VI Case Study

VI-A System Description

In order to analyze the performance of the stochastic UC model presented in Section V, we investigate a modified version of the IEEE RTS-96 power system from [20] depicted in Fig. 5, with 48 buses comprising areas 1 and 2 of the original system. Table II shows the relevant parameters of different thermal plant types. The studied system includes 20 generators and 16 wind farms. We assume that six wind farms are providing virtual inertia; four of them via VSM control, and the remaining two through equivalent droop regulation. The UC is ran for two days without frequency constraints in order to initialize the system prior to introducing the frequency constraints on days 3, 4 and 5. This is done to ensure the impact of start-up costs are well distributed and not concentrated on one day. Therefore, the total simulation horizon is five days ($T=120\,\mathrm{h}$) whereas the UC schedule is optimized separately for each day, with the last hour of each day used as an input for the next.

As the set of possible contingencies we consider the failure of synchronous generators $i=\{1,6,8,10\}$. These generators are of various capacities, ranging from the smallest to the largest unit in the system. The hour 19 of day 3 (i.e., $t=67\,\mathrm{h}$) is selected to be the time instance of a possible contingency, as this is the hour with high wind penetration and low demand. Ten wind power scenarios are considered, which brings the total number of scenarios to 50. The optimization problem is formulated in Python and uses the Gurobi solver with default parameterization.

VI-B Results

In this section, the simulation results from the stochastic UC are presented. Fig. 6 showcases the load and wind power profiles as well as the aggregate dispatch of synchronous generation for two UC runs: (i) without frequency constraints; and (ii) with frequency constraints. Furthermore, Table III indicates the difference in the total number of generators committed between the two runs. Both Fig. 6 and Table III suggest that, although the amount of committed generators increases significantly, the total SG production is only slightly changed. This is justified by the fact that the extra generators are solely committed for the purpose of providing inertia, and are thus operating at their technical minimum. The production surplus arising from the additionally committed units is compensated by wind curtailment and other generators reducing their power output.

Moreover, the evolution of aggregate system inertia over the course of the whole scheduling horizon is depicted in Fig. 7. A noticeable step change in total system inertia at hour $67$ reflects the violation of frequency constraints under contingency, which subsequently triggers a dispatch of auxiliary synchronous generators. While the inertia levels do not differentiate between the two scenarios during the first two days, on the days following the potential outage some carryover impacts can be observed. This is a consequence of the commitment schedule being radically changed at hour $t=67$, thus affecting the UC schedules in the following days.

Some insightful conclusions can be drawn from Fig. 8, where the difference between the actual values of the frequency metrics and the respective ENTSO-E thresholds are depicted. For this purpose, we define a constraint gap $\eta$ as a measure of the relative constraint distance to its limit, e.g., $\eta_{\mathrm{nadir}}=\Delta f_{\mathrm{max}}/\Delta f_{\mathrm{lim}}-1$. After the completion of unit commitment, the constraints are re-evaluated using the obtained $F_{g}$, $R_{g}$ and $M$ values in order to determine which frequency criteria becomes binding at the instance of the fault. A negative constraint gap corresponds to the non-binding constraint, i.e., the specific frequency criteria being met. It should be noted that the positive values of $\eta$ for $t<67\,\mathrm{h}$ indicate that the frequency threshold would be violated if a fault occurs. No action is required however, considering that in this case study we assume that the contingency can only occur in hour $t=67$. The Fig. 8 indicates that without explicit frequency metric constraints all of these constraints would be violated. Moreover, it can be observed that when including the frequency constraints, the RoCoF constraint is closest to its limit - corresponding to the smallest constraint gap - and thus binding. The constraint gap difference between the two scenarios at hour $t=67$ clearly highlights the importance of including the frequency constraints in UC in order to avoid large frequency excursions and undesired triggering of protection and UFLS schemes. The same observations are also reflected in Fig. 9 through time-domain frequency response of the system. Understandably, the values of RoCoF, nadir and steady-state deviation are reduced compared to the scenario without frequency constraints, such that all of the ENTSO-E criteria are fulfilled.

Finally, we investigate the economic impact of including the frequency constraints into the stochastic UC model. The breakdown of operational costs for day 3 is presented in Table IV. The addition of frequency constraints leads to a $5\,\%$ increase in total expected system costs and a significant increase in start-up costs by $185\,\%$. This is due to six extra generators being turned on for providing inertia at the period of a potential generation failure, as shown in Table III. A large increase is also seen in reserve scheduling costs, as the reserves are now not only scheduled to cover wind power uncertainty but also for possible contingencies. The change in cost is of course highly dependent on the specific system and the considered contingencies.

VII Conclusion

This paper includes frequency constraints in the UC problem of a system with large wind power penetration in order to investigate the impact of frequency dynamics on unit scheduling. By employing the analytic expressions for post-contingency frequency response of a multi-machine system, we define a set of constraints reflecting the frequency nadir, RoCoF and quasi steady-state deviation. The highly non-linear frequency nadir constraint is linearized using two approaches: (i) a PWL technique adapted from the literature; and (ii) a proposed simple and efficient method for extracting bounds on decision variables of interest, which is shown to be computationally superior to PWL. Using the latter approach, the stochastic UC problem is formulated as an MILP, with an objective of minimizing the expected system costs against wind power production and generation outage uncertainties.

Our results show that the inclusion of frequency constraints in the UC model significantly affects the dispatch of synchronous generators and consequently the expected system costs. Indeed, during anticipated critical events such as the loss of generation, additional synchronous machines are needed for providing sufficient inertia and damping in the process of frequency containment. Such actions lead to a drastic increase in the UC costs, especially start-up and reserve scheduling, which poses a new challenge as the operator must find a way to remunerate the units committed for the sole purpose of frequency regulation. This is an exciting avenue for future work.

Appendix A

For clarity, we visually illustrate in Fig. 10 how the PWL optimization problem from (6) is solved. In this example there are four evaluation points set at -7.5, -2.5, 2.5 and 7.5, and the respective function is approximated with four segments. At each evaluation point the model identifies the segment closest to the original curve, and subsequently aims to minimize the overall shaded area. The proposed technique can then be expanded and employed on a function of three variables, as we have done with the expression for frequency nadir.

References section