Multi-stage Robust Unit Commitment with Discrete Load Shedding Based on Partially Affine Policy and Two-stage Reformulation

Zhongjie Guo; Jiayu Bai; Wei Wei; Haifeng Qiu; Weihao Hu

网刊加载中。。。

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

确定继续浏览么?

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

Multi-stage Robust Unit Commitment with Discrete Load Shedding Based on Partially Affine Policy and Two-stage Reformulation PDF

- ORCID：
Zhongjie Guo ¹ (Member, IEEE)
✉
- ORCID：
Jiayu Bai ² (Member, IEEE)
✉
- ORCID：
Wei Wei ³ (Senior Member, IEEE)
✉
- ORCID：
Haifeng Qiu ⁴ (Member, IEEE)
✉
- ORCID：
Weihao Hu ¹ (Senior Member, IEEE)
✉

1. School of Mechanical and Electrical Engineering, University of Electronic Science and Technology of China, Chengdu 610054, China； 2. State Grid Sichuan Electric Power Research Institute, Chengdu 610041, China； 3. Department of Electrical Engineering, Tsinghua University, Beijing 100084, China； 4. School of Electrical and Electronic Engineering, Nanyang Technological University, Singapore 639798, Singapore

Updated：2025-03-26

DOI：10.35833/MPCE.2024.000202

OUTLINE

Abstract

This paper studies the problem of multi-stage robust unit commitment with discrete load shedding. In the day-ahead phase, the on-off status of thermal units is scheduled. During each period of real-time dispatch, the output of thermal units and the action of load shedding are determined, and the discrete choice of load shedding corresponds to the practice of tripping substation outlets. The entire decision-making process is formulated as a multi-stage adaptive robust optimization problem with mixed-integer recourse, whose solution takes three steps. First, we propose and apply partially affine policy, which is optimized ahead of the day and restricts intertemporal dispatch variables as affine functions of previous uncertainty realizations, leaving remaining continuous and binary dispatch variables to be optimized in real time. Second, we demonstrate that the resulting model with partially affine policy can be reformulated as a two-stage robust optimization problem with mixed-integer recourse. Third, we modify the standard nested column-and-constraint generation algorithm to accelerate the inner loops by warm start. The modified algorithm solves the two-stage problem more efficiently. Case studies on the IEEE 118-bus system verify that the proposed partially affine policy outperforms conventional affine policy in terms of optimality and robustness; the modified nested column-and-constraint generation algorithm significantly reduces the total computation time; and the proposed method balances well optimality and efficiency compared with state-of-the-art methods.

Keywords

Unit commitment; robust optimization; uncertainty; affine policy; load shedding

I. Introduction

RENEWABLE energy generation introduces unprecedented randomness to power system operations, so the grid operator should make decisions considering how to handle uncertainty. Unit commitment (UC) is one of the most important decision-making problems. Using robust optimization theory, two-stage robust UC has been widely investigated [

1]-[3], which divides the variables into two groups. One group includes commitment variables, which are called here-and-now decisions and determined before uncertainty observation. The other group includes all the dispatch variables, which are called wait-and-see recourse variables and determined after uncertainty observation.

A shortcoming of two-stage formulation is pointed out via a counterexample in [

4] that the non-anticipativity of dispatch decisions is violated. Non-anticipativity is a crucial issue in decision-making problems with uncertainty, which refers to that the decisions in any period (stage) cannot rely on the uncertainty realizations in the future [5]. In other words, the underlying decision rule of

x_{t}

should be

x_{t} = f_{t} (ξ_{t})

where

ξ_{t}

is the uncertainty realizations up to period

t

. The dispatch sequence behind the two-stage formulation is

\{(ξ_{1}, ξ_{2}, \dots, ξ_{T}), (x_{1}, x_{2}, \dots, x_{T})\}

, which does not respect the non-anticipativity.

To address this issue, multi-stage robust UC (MRUC) has been studied, which formulates the real-time dispatch and uncertainty realization as two intertwined processes, i.e., $\{ξ_{1}, x_{1}, ξ_{2}, x_{2}, \dots, ξ_{T}, x_{T}\}$ , which aligns with operation practice [

6]. The MRUC problem has a nested structure and thus is intractable. Most existing research works consider continuous dispatch variables with three mainstream solution methods.

1) Policy/rule-based methods [

7], [8]. These methods stipulate the dispatch variable

x_{t}

as a specific function (typically affine function) of previous uncertainty realizations. The affine policy is

x_{t} = a_{t} ξ_{t} + b_{t},

where

a_{t}

and

b_{t}

are the coefficient variables that are determined in the day-ahead phase. Such methods naturally satisfy non-anticipativity because the rules have specified the dependence on previous information only. Using the affine policy, the resulting model is a static robust optimization problem and can be solved by its robust counterpart or other techniques. However, the pre-selected affine structure may jeopardize the optimality.

2) Region-based methods, also named as implicit policy method [

9]-[11]. These methods aim to find time-independent feasible regions/bounds for intertemporal dispatch variables, which are established ahead of the day by considering some vertices of uncertainty set and extreme ramping events between two consecutive periods. With these regions, the real-time dispatch resorts to optimal power flow problems that are decoupled over time.

3) Dynamic programming (DP)-based methods [

12], [13]. Following the convention of DP, the optimization problem in each stage minimizes the sum of instant cost and future cost. The future cost is estimated by cost-to-go functions, which are trained ahead of the day. Due to the curses of dimensionality, DP-based methods suffer from high computation burden; besides, robust feasibility is a tough issue for DP-based methods, so they usually assume a relatively complete recourse, which is impossible in practice.

The above methods handle MRUC problem with continuous dispatch variables in their own ways. In practice, some real-time dispatch actions are discrete, such as load shedding. Existing research works treat load shedding as a continuous variable, which should be discrete [

14]. The reason is that load shedding is usually realized by tripping substation outlets [15]. From the perspective of transmission system, the load cannot be continously adjusted. As a result, we need to consider MRUC with mixed-integer dispatch variables.

With integer recourse, the problem becomes much more complicated. If we revisit the aforementioned three methods, the region-based methods do not work because they rely on the model linearity and convexity. The DP-based method is extended to consider integer variables in [

16], but the computation burden remains a trouble due to curses of dimensionality and higher complexity. Besides, to use decision rules, one task is to tackle binary recourse, which requires piecewise constant policy. Finite adaptability technique makes sense in two-stage problems [17], [18], but no successful attempt is found in MRUC applications. The state-of-the-art rule-based methods are reported in [19]-[21]. The method in [19] and [20] uses a pre-defined partitioning function to formulate the piecewise constant rules, and the final problem is a scalable mixed-integer linear program (MILP). However, there is no specific method to design the optimal partitioning function. Instead, the method in [21] optimizes the partitioning function ahead of the day, but the scalability is jeopardized.

Therefore, this paper focuses on the MRUC problem with discrete load shedding and investigates how to handle binary recourse in a more tractable manner. The contributions are twofold below.

1) Problem modeling: we formulate the UC with discrete load shedding as a multi-stage adaptive robust optimization problem with mixed-integer recourse. The discrete choice of load shedding corresponds to the practice of tripping substation outlets. In the multi-stage decision-making sequence, the on-off status of thermal units is scheduled ahead of the day. During each period of real-time dispatch, the output of thermal units and the action of load shedding are determined.

2) Solution methodology: we propose a systematic solution method to MRUC problem. The first step is to apply partially affine policy to restrict intertemporal dispatch variables as affine functions of previous uncertainty realizations. Non-state and binary dispatch variables are left to be optimized in real time. The second step is to derive a two-stage robust optimization problem with partially affine policy. The third step is to implement a modified nested column-and-constraint generation (M-nCCG) algorithm, which refines the standard nested column-and-constraint generation (S-nCCG) algorithm by warming up the inner loops.

The rest of this paper is organized as follows. The mathematical formulation is introduced in Section II. The multi-stage adaptive robust optimization model is presented in Section III, whose solution methodology is proposed in Section IV. Case studies are provided in Section V, and conclusions are drawn in Section VI.

II. Mathematical Formulation

A. Thermal Units

1) Coal-fired Units

b_{t}^{c}, u_{t}^{c}, v_{t}^{c} \in \{0,1\} \forall c \in C, t \in T

(1a)

- b_{t - 1}^{c} + b_{t}^{c} - b_{τ}^{c} \leq 0 \forall c \in C, t \in T, τ = (t + 1) : m i n \{β_{c}^{U} + t - 1, T\}

(1b)

b_{t - 1}^{c} - b_{t}^{c} + b_{τ}^{c} \leq 1 \forall c \in C, t \in T, τ = (t + 1) : m i n \{β_{c}^{D} + t - 1, T\}

(1c)

b_{t - 1}^{c} - b_{t}^{c} + u_{t}^{c} - v_{t}^{c} = 0 \forall c \in C, t \in T

(1d)

P_{t}^{c} - P_{t - 1}^{c} \leq R_{c}^{+} (1 + b_{t - 1}^{c} - b_{t}^{c}) + (2 - b_{t}^{c} - b_{t - 1}^{c}) L^{c} \forall c \in C, t \in T

(1e)

P_{t - 1}^{c} - P_{t}^{c} \leq R_{c}^{-} (1 - b_{t - 1}^{c} + b_{t}^{c}) + (2 - b_{t}^{c} - b_{t - 1}^{c}) U^{c} \forall c \in C, t \in T

(1f)

b_{t}^{c} P_{m i n}^{c} \leq P_{t}^{c} \leq b_{t}^{c} P_{m a x}^{c} \forall c \in C, t \in T

(1g)

where $c$ and $C$ are the index and set of coal-fired units, respectively; $t$ and $T = \{1,2, \dots, T\}$ are the index and set of periods, respectively; $b_{t}^{c}$ , $u_{t}^{c}$ , $v_{t}^{c}$ , and $P_{t}^{c}$ are the decision variables, and if $b_{t}^{c} = 1$ ( $b_{t}^{c} = 0$ ), unit $c$ is on (off) during period $t$ , and $u_{t}^{c} = 1$ ( $v_{t}^{c} = 1$ ) means unit $c$ is turned on (off) during period $t$ ; $β_{c}^{U}$ and $β_{c}^{D}$ are the minimum-up time and minimum-down time, respectively; $R_{c}^{+}$ and $R_{c}^{-}$ are the ramp-up and ramp-down rates, respectively, meanwhile, if unit $c$ starts up during period $t$ , its output $P_{t}^{c}$ is no higher than the start-up rate $L^{c}$ , and if it shuts down during period $t$ , $P_{t - 1}^{c}$ is no higher than the shut-down rate $U^{c}$ ; and $P_{m i n}^{c}$ and $P_{m a x}^{c}$ are the lower and upper bounds on $P_{t}^{c}$ , respectively. Formulas (1b) and (1c) are the minimum-up and minimum-down time constraints, respectively; if unit $c$ starts up during period $t$ , i.e., $b_{t - 1}^{c} = 0$ and $b_{t}^{c} = 1$ , (1b) keeps this unit on during the next $β_{c}^{U} - 1$ periods, which is similar to the minimum-down time. According to (1d), we have $u_{t}^{c} = 1$ if $b_{t - 1}^{c} = 0$ and $b_{t}^{c} = 1$ , and $v_{t}^{c} = 1$ if $b_{t - 1}^{c} = 1$ and $b_{t}^{c} = 0$ . If $b_{t - 1}^{c} = b_{t}^{c}$ , both $u_{t}^{c}$ and $v_{t}^{c}$ will be zero since the start-up and shut-down costs are to be minimized. Formulas (1e) and (1f) are the improved ramping constraints. Finally, (1g) imposes upper and lower bounds on $P_{t}^{c}$ .

2) Gas-fired Units

Compared with coal-fired units, the gas-fired ones can reach the maximum output within a few minutes. Therefore, in an hourly model, the ramping limits can be discarded [

22]. The operating model of gas-fired units is given below.

b_{t}^{g}, u_{t}^{g}, v_{t}^{g} \in \{0,1\} \forall g \in G, t \in T

(2a)

- b_{t - 1}^{g} + b_{t}^{g} - b_{τ}^{g} \leq 0 \forall g \in G, t \in T, τ = (t + 1) : m i n \{β_{g}^{U} + t - 1, T\}

(2b)

b_{t - 1}^{g} - b_{t}^{g} + b_{τ}^{g} \leq 1 \forall g \in G, t \in T, τ = (t + 1) : m i n \{β_{g}^{D} + t - 1, T\}

(2c)

b_{t - 1}^{g} - b_{t}^{g} + u_{t}^{g} - v_{t}^{g} = 0 \forall g \in G, t \in T

(2d)

b_{t}^{g} P_{m i n}^{g} \leq P_{t}^{g} \leq b_{t}^{g} P_{m a x}^{g} \forall g \in G, t \in T

(2e)

where $g$ and $G$ are the index and set of gas-fired units, respectively. We do not interpret model (2) in detail, because it can be easily understood by analogy with model (1). The main difference is that the gas-fired unit is free from ramping limits. The minimum-up and minimum-down constraints remain, but $β_{g}^{U}$ and $β_{g}^{D}$ may be smaller. In some research works, the ramping rate of gas-fired units is considered [

23], [24]. Nevertheless, we can always divide all thermal units into ramping-constrained ones and ramping-free ones in an hourly UC problem.

B. Uncertain Wind Power

Let $w$ and $W$ be the index and set of wind farms, respectively. Under time-varying weather conditions, the available wind power $ξ_{t}^{w}$ is random. The utilized wind power is denoted by $P_{t}^{w}$ , and the rest $Δ P_{t}^{w}$ is curtailed. The operating model of wind generation is given below.

\{\begin{array}{l} 0 \leq P_{t}^{w} \leq ξ_{t}^{w} \\ P_{t}^{w} + Δ P_{t}^{w} = ξ_{t}^{w} \end{array} \forall w \in W, \forall t \in T

(3)

To handle the randomness of wind power, we use the following box-type (polyhedral) uncertainty set:

\begin{array}{l} Ξ_{t}^{w} = \{ξ_{t}^{w} | ξ_{t}^{w, 0} - δ_{t}^{w} \leq ξ_{t}^{w} \leq ξ_{t}^{w, 0} + δ_{t}^{w}, 0 \leq ξ_{t}^{w} \leq P_{m a x}^{w}\} \\ \forall w \in W, \forall t \in T \end{array}

(4)

where $ξ_{t}^{w, 0}$ is the wind power forecast; $δ_{t}^{w}$ is the maximum forecast error; and $P_{m a x}^{w}$ is the installed capacity of wind generation. If necessary, spatial budget can be incorporated into the uncertainty set. The resulting uncertainty set remains a polyhedron, which can be handled by the method to be proposed. However, the temporal budget is usually not considered in the MRUC problem [

10], because it contains global information that will violate the non-anticipativity of decisions.

C. Load

Let $l$ and $L$ be the index and set of loads, respectively. We always hope that load $P_{t}^{l}$ can be fully satisfied, but load shedding $Δ P_{t}^{l}$ is sometimes inevitable. In the real world, load shedding is usually realized by tripping substation outlets [

15]. As in Fig. 1,

P_{t}^{l}

is distributed in

N_{l}

outlets. In practice, the load on each outlet is measured. Without loss of generality, we assume

P_{t}^{l, 1} = P_{t}^{l, 2} = \dots = P_{t}^{l, N_{l}} = P_{t}^{l} / N_{l}

Fig. 1 Load shedding by tripping substation outlet.

We assume the former $N_{l}^{*}$ outlets can be tripped. Outlets $(N_{l}^{*} + 1) : N_{l}$ connecting to important loads, such as hospital, should not be tripped. Therefore, using binary variables $b_{t}^{l, 1 : N_{l}^{*}}$ , the load model is:

b_{t}^{l, n} \in \{0,1\} n = 1 : N_{l}^{*}, \forall l \in L, \forall t \in T

(5a)

Δ P_{t}^{l} = \frac{P_{t}^{l}}{N_{l}} \sum_{n = 1}^{N_{l}^{*}} b_{t}^{l, n} \forall l \in L, \forall t \in T

(5b)

Like much UC research, load $P_{t}^{l}$ is given as a parameter [

25]. Nevertheless, the proposed method can handle load uncertainty in the same way as wind uncertainty.

Two remarks are given below.

1) Load shedding is usually a small probability event. However, making UC decisions with load shedding is meaningful in power systems with a high penetration of renewable generation. To respond the fluctuation of renewable power, strategic load shedding can be regarded as an emergency measure, ensuring the system to continue running with the minimum negative effects.

2) The substation in Fig. 1 interfaces the transmission network with downstream distribution networks. In other words, the $n$ ^th outlet is the inlet of the $n$ ^th distribution substation. If necessary, the outlets of distribution substations can be considered in model (5). As a result, $N_{l}$ will be the total number of distribution outlets. Therefore, the proposed method can unify transmission and distribution systems, optimizing the power balance in a wider spatial range.

D. Network Constraints

The renowned direct-current power flow is employed in the form of power transfer distribution factor (PTDF):

\begin{matrix} \sum_{\forall c \in C} P_{t}^{c} + \sum_{\forall g \in G} P_{t}^{g} + \sum_{\forall w \in W} P_{t}^{w} = \sum_{\forall l \in L} (P_{t}^{l} - Δ P_{t}^{l}) \forall t \in T \end{matrix}

(6a)

\begin{array}{l} - P_{m a x}^{f} \leq \sum_{\forall c \in C} π_{c f} P_{t}^{c} + \sum_{\forall g \in G} π_{g f} P_{t}^{g} + \sum_{\forall w \in W} π_{w f} P_{t}^{w} - \\ \sum_{\forall l \in L} π_{l f} (P_{t}^{l} - Δ P_{t}^{l}) \leq P_{m a x}^{f} \end{array}

(6b)

where $f$ and $F$ are the index and set of transmission lines, respectively; $π_{(\cdot) f}$ is the PTDF from a certain facility to line $f$ ; and $P_{m a x}^{f}$ is the transmission capacity. Equation (6a) is the system-wise power balance condition. Formula (6b) bounds the active power flow on each line $f$ .

E. Cost Function

The costs in UC include the start-up/shun-down cost and fuel cost of thermal units, as well as the penalty of wind energy curtailment and load shedding. Hence, the cost function is expressed as:

\begin{array}{l} \sum_{\forall t \in T} [\sum_{\forall c \in C} (C_{c}^{S U} u_{t}^{c} + C_{c}^{S D} v_{t}^{c} + C_{c}^{f i x e d} b_{t}^{c}) \sum_{\forall g \in G} (C_{g}^{S U} u_{t}^{g} + C_{g}^{S D} v_{t}^{g} + \\ \overset{}{C_{g}^{f i x e d} b_{t}^{g}})] + \sum_{\forall t \in T} h (\sum_{\forall c \in C} C_{c}^{F u e l} P_{t}^{c} + \sum_{\forall g \in G} C_{g}^{F u e l} P_{t}^{g} + \\ \sum_{\forall l \in L} C_{l}^{S h e d} Δ P_{t}^{l} + \sum_{\forall w \in W} C_{w}^{C u r t} Δ P_{t}^{w}) \end{array}

(7)

where $C_{c}^{S U}$ , $C_{g}^{S U}$ , $C_{c}^{S D}$ , and $C_{g}^{S D}$ are the costs of start-up $a n d$ shut-down, respectively; $C_{c}^{f i x e d}$ and $C_{g}^{f i x e d}$ are the fixed costs of running for one hour; $C_{c}^{F u e l}$ and $C_{g}^{F u e l}$ are the fuel cost coefficients; $C_{l}^{S h e d}$ and $C_{w}^{C u r t}$ are the penalty coefficients of load shedding and wind energy curtailment, respectively, and $C_{l}^{S h e d}$ is much larger; and constant $h$ is 1 hour. One can also use quardratic fuel functions, which can be approximated by linear pieces [

26].

III. MRUC with Discrete Load Shedding

A. Notations

For brevity, we define some notations. All decision variables are encapsulated into commitment vector $φ$ , continuous intertemporal (state) vector $x_{t}$ , continuous non-state vector $y_{t}$ , binary non-state vector $z_{t}$ , and random vector $ξ_{t}$ .

φ = {[φ_{1}, φ_{2}, . . ., φ_{|T|}]}^{T}

(8a)

φ_{t a} = [b_{t}^{c_{1}}, u_{t}^{c_{1}}, v_{t}^{c_{1}}, b_{t}^{c_{2}}, u_{t}^{c_{2}}, v_{t}^{c_{2}}, . . ., b_{t}^{c_{|C|}}, u_{t}^{c_{|C|}}, v_{t}^{c_{|C|}}]

(8b)

φ_{t b} = [b_{t}^{g_{1}}, u_{t}^{g_{1}}, v_{t}^{g_{1}}, b_{t}^{g_{2}}, u_{t}^{g_{2}}, v_{t}^{g_{2}}, . . ., b_{t}^{g_{|C|}}, u_{t}^{g_{|C|}}, v_{t}^{g_{|C|}}]

(8c)

φ_{t} = {[φ_{t a}, φ_{t a}]}^{}

(8d)

x_{t} = {[P_{t}^{c_{1}}, P_{t}^{c_{2}}, . . ., P_{t}^{c_{|C|}}]}^{T} \forall t \in T

(8e)

\{\begin{array}{l} P_{}^{g} = {[P_{t}^{g_{1}}, P_{t}^{g_{2}}, . . ., P_{t}^{g_{|G|}}]}^{T} \\ P_{}^{w} = {[P_{t}^{w_{1}}, P_{t}^{w_{2}}, . . ., P_{t}^{w_{|W|}}]}^{T} \\ Δ P_{}^{w} = {[Δ P_{t}^{w_{1}}, Δ P_{t}^{w_{1}}, . . ., Δ P_{t}^{w_{|W|}}]}^{T} \\ Δ P_{}^{l} = {[Δ P_{t}^{l_{1}}, Δ P_{t}^{l_{2}}, . . ., Δ P_{t}^{l_{|L|}}]}^{T} \\ y_{t} = {[P_{}^{g}, P_{}^{w}, Δ P_{}^{w}, Δ P_{}^{l}]}^{T} \end{array} \forall t \in T

(8f)

b_{t}^{l_{1}} = [b_{t}^{l_{1}, 1}, b_{t}^{l_{1}, 2}, . . ., b_{t}^{l_{1}, N_{1}^{*}}]

(8g)

b_{t}^{l_{2}} = [b_{t}^{l_{2}, 1}, b_{t}^{l_{2}, 2}, . . ., b_{t}^{l_{2}, N_{2}^{*}}]

(8h)

b_{t}^{l_{|L|}} = [b_{t}^{l_{|L|}, 1}, b_{t}^{l_{|L|}, 2}, . . ., b_{t}^{l_{|L|}, N_{|L|}^{*}}]

(8i)

z_{t} = {[b_{t}^{l_{1}}, b_{t}^{l_{2}}, . . ., b_{t}^{l_{|L|}}]}^{T}

(8j)

ξ_{t} = {[ξ_{t}^{1}, ξ_{t}^{2}, . . ., ξ_{t}^{| 𝕎 |}]}^{T} \forall t \in T

(8k)

Then, we define the feasible sets as:

Φ : = \{φ |(1 a) - (1 d), (2 a) - (2 d)\}

(9a)

Ξ_{t} = \{Ξ_{t}^{1}, Ξ_{t}^{2}, . . ., Ξ_{t}^{| 𝕎 |}\} \forall t \in T

(9b)

Ω_{t} (φ, ξ_{t}, x_{t - 1}) : = \{(x_{t}, y_{t}, z_{t})| (1 e), (1 f), (1 g), (2 e), (3), (5), (6)\} \forall t \in T

(9c)

All the constraints in $Ω_{t}$ are linear, so we can rewrite it as:

Ω_{t} (φ, ξ_{t}, x_{t - 1}) : = \{(x_{t}, y_{t}, z_{t}) | A_{t} x_{t} + B_{t} y_{t} + C_{t} z_{t} \leq b_{t} + D_{t} φ + E_{t} ξ_{t} + F_{t} x_{t - 1}\}

(9d)

where $A_{t}$ , $B_{t}$ , $C_{t}$ , $D_{t}$ , $E_{t}$ , $F_{t}$ , and $b_{t}$ are the constant matrices extracted from (1e), (1f), (1g), (2e), (3), (5), and (6), respectively.

Besides, the cost function in (7) is expressed as:

k^{T} φ + \sum_{t = 1}^{T} (r^{T} x_{t} + s^{T} y_{t})

where $k^{T} φ$ is the UC cost; $\sum_{t = 1}^{T} (r^{T} x_{t} + s^{T} y_{t})$ is the total dispatch cost; and $k$ , $r$ , and $s$ are the coefficient vectors.

B. MRUC

When randomness is absent, UC is a deterministic optimization problem, which minimizes the cost in (7) subject to the constraints in (1)-(3), (5), and (6). However, the wind power is uncertain. The entire decision-making sequence is:

\underset{D a y - a h e a d}{\underset{︸}{φ}}, \underset{R e a l - t i m e}{\underset{︸}{ξ_{1}, (x_{1}, y_{1}, z_{1}), ξ_{2}, (x_{2}, y_{2}, z_{2}), \dots, ξ_{T}, (x_{T}, y_{T}, z_{T})}}

(10)

In the day-ahead phase, UC $φ$ is determined. During each period $t$ of real-time dispatch, decisions $(x_{t}, y_{t}, z_{t})$ are made after the wind power $ξ_{t}$ is observed. Such a sequence is non-anticipated, i.e., dispatch decisions during any period $t$ do not rely on the uncertainty realizations in the future.

To describe (10) and manage the uncertainty in a robust manner, we formulate the MRUC problem as:

\underset{φ \in Φ}{m i n} k^{T} φ + (\underset{ξ_{1} \in Ξ_{1}}{m a x} \underset{(x_{1}, y_{1}, z_{1}) \in Ω_{1} (φ, ξ_{1}, x_{0})}{m i n} r^{T} x_{1} + s^{T} y_{1}) + (\underset{ξ_{2} \in Ξ_{2}}{m a x} \underset{(x_{2}, y_{2}, z_{2}) \in Ω_{2} (φ, ξ_{2}, x_{1})}{m i n} r^{T} x_{2} + s^{T} y_{2}) + \dots + (\underset{ξ_{T} \in Ξ_{T}}{m a x} \underset{(x_{T}, y_{T}, z_{T}) \in Ω_{T} (φ, ξ_{T}, x_{T - 1})}{m i n} r^{T} x_{T} + s^{T} y_{T})

(11)

IV. Solution Methodology

The MRUC problem (11) is intractable for two reasons. One is the nested structure coupled by intertemporal variables, and the other is the non-convexity and non-continuity caused by binary recourse. To solve (11), we take three steps: first, propose a partially affine policy and apply it to (11); second, establish an equivalent two-stage robust optimization problem with mixed-integer recourse; third, solve the equivalent problem by an M-nCCG algorithm.

A. Partially Affine Policy

Fully affine policy is proposed in [

27], restricting continuous dispatch decisions during period

t

to be affine functions of previous uncertainty realizations. Define

ξ_{[t]} = [ξ_{1}, ξ_{2}, \dots, ξ_{t}]

, and the fully affine policy is:

[x_{t}; y_{t}] = U_{t} ξ_{[t]} + V_{t} \forall t \in T

(12)

where $U_{t}$ and $V_{t}$ are the decision matrices that will be optimized together with $φ$ in the day-ahead phase. However, there are two main limitations:

1) Fully affine policy is criticized for suboptimality. Especially in the presence of gas-fired units which take non-state real-time actions, the affine relation jeopardizes the real-time flexibility.

2) Affine policy cannot handle binary recourse. The latest method in [

28] establishes piecewise constant polices for binary recourse

z_{t}

. However, this method is unfriendly for its complicated partitioning and lifting techniques, as well as the priori and arbitrary choice of a non-linear function.

To address the first limitation, we introduce partially affine policy to unleash non-state variables and impose affine relation only on intertemporal dispatch variables:

x_{t} = U_{t} ξ_{[t]} + V_{t} \forall t \in T

(13)

During each period $t$ of real-time dispatch, $x_{t}$ is determined based on (13) after $ξ_{[t]}$ is observed. How to determine $y_{t}$ and $z_{t}$ will be discussed later in Section IV-D.

B. Two-stage Equivalence

We define some new notations marked by overlines as: $\bar{φ} : = {[φ^{T}, U_{1}, V_{1}, U_{2}, V_{2}, . . ., U_{|T|}, V_{|T|}]}^{T}$ , $\bar{Φ} = \{\bar{φ} | φ \in Φ\}$ and ${\bar{Ω}}_{t} (\bar{φ}, ξ_{[t]}) : =$ $\{(y_{t}, z_{t}) | B_{t} y_{t} + C_{t} z_{t} \leq \begin{matrix} {\bar{b}}_{t} \end{matrix}\}$ , where $\begin{matrix} {\bar{b}}_{t} \end{matrix} = b_{t} + D_{t} φ + E_{t} ξ_{t} + F_{t} (U_{t - 1} ξ_{[t - 1]} + V_{t - 1}) - A_{t} (U_{t} ξ_{[t]} + V_{t})$ .

By applying partially affine policy (13) to problem (11), we derive:

\begin{array}{l} \underset{\bar{φ} \in \bar{Φ}}{m i n} (k^{T} φ + \sum_{t = 1}^{T} r^{T} V_{t}) + (\underset{ξ_{1} \in Ξ_{1}}{m a x} \underset{(y_{1}, z_{1}) \in {\bar{Ω}}_{1} (\bar{φ}, ξ_{[1]})}{m i n} r^{T} U_{1} ξ_{[1]} + s^{T} y_{1}) + \\ (\underset{ξ_{2} \in Ξ_{2}}{m a x} \underset{(y_{2}, z_{2}) \in {\bar{Ω}}_{2} (\bar{φ}, ξ_{[2]})}{m i n} r^{T} U_{2} ξ_{[2]} + s^{T} y_{2}) + \dots + (\underset{ξ_{T} \in Ξ_{T}}{m a x} \underset{(y_{T}, z_{T}) \in {\bar{Ω}}_{T} (\bar{φ}, ξ_{[T]})}{m i n} r^{T} U_{T} ξ_{[T]} + s^{T} y_{T}) \end{array}

(14)

Proposition 1 Problem (14) is equivalent to the following two-stage robust optimization problem with mixed-integer recourse:

\underset{\bar{φ} \in \bar{Φ}}{m i n} (k^{T} φ + \sum_{t = 1}^{T} r^{T} V_{t}) + \underset{ξ_{t} \in Ξ_{t}, \forall t \in T}{m a x} \underset{(y_{t}, z_{t}) \in {\bar{Ω}}_{t} (\bar{φ}, ξ_{[t]}), \forall t \in T}{m i n} \sum_{t = 1}^{T} (r^{T} U_{t} ξ_{[t]} + s^{T} y_{t})

(15)

Proof Below is the dispatch during periods $T - 1$ and $T$ .

\begin{array}{l} (\underset{ξ_{T - 1} \in Ξ_{T - 1}}{m a x} \underset{(y_{T - 1}, z_{T - 1}) \in {\bar{Ω}}_{T - 1} (\bar{φ}, ξ_{[T - 1]})}{m i n} r^{T} U_{T - 1} ξ_{[T - 1]} + s^{T} y_{T - 1}) + \\ (\underset{ξ_{T} \in Ξ_{T}}{m a x} \underset{(y_{T}, z_{T}) \in {\bar{Ω}}_{T} (\bar{φ}, ξ_{[T]})}{m i n} r^{T} U_{T} ξ_{[T]} + s^{T} y_{T}) \end{array}

(16)

Notice that the minimization over $(y_{T - 1}, z_{T - 1})$ is independent from the maximization over $ξ_{T}$ and the minimization over $(y_{T}, z_{T})$ , implying that the optimal solution of (16) remains the same if we exchange the inner minimum and maximum operators [

29]. Therefore, (16) is equivalent to (17a) and by aggregation further equivalent to (17b).

\begin{array}{l} (\underset{ξ_{T - 1} \in Ξ_{T - 1}}{m a x} \underset{ξ_{T} \in Ξ_{T}}{m a x} \underset{(y_{T - 1}, z_{T - 1}) \in {\bar{Ω}}_{T - 1} (\bar{φ}, ξ_{[T - 1]})}{m i n} r^{T} U_{T - 1} ξ_{[T - 1]} + s^{T} y_{T - 1}) + \\ (\underset{(y_{T}, z_{T}) \in {\bar{Ω}}_{T} (\bar{φ}, ξ_{[T]})}{m i n} r^{T} U_{T} ξ_{[T]} + s^{T} y_{T}) \end{array}

(17a)

\underset{\begin{matrix} ξ_{T - 1} \in Ξ_{T - 1}, \\ ξ_{T} \in Ξ_{T} \end{matrix}}{m a x} \underset{\begin{matrix} (y_{T - 1}, z_{T - 1}) \in {\bar{Ω}}_{T - 1} (\bar{φ}, ξ_{[T - 1]}), \\ (y_{T}, z_{T}) \in {\bar{Ω}}_{T} (\bar{φ}, ξ_{[T]}) \end{matrix}}{m i n} \sum_{t = T - 1}^{T} (r^{T} U_{t} ξ_{[t]} + s^{T} y_{t})

(17b)

In a similar way, combining the dispatch during period $T - 2$ with (17) will produce a maximum-minimum equivalence involving periods $T - 2$ , $T - 1$ , and $T$ . By a backward induction to the first period, we can derive (15), which completes the proof.

C. M-nCCG Algorithm

According to Proposition 1, solving MRUC problem (14) is equivalent to solving (15), which is a two-stage robust optimization problem with mixed-integer recourse. The mainstream solution method is the nested column-and-constraint generation (CCG) algorithm proposed in [

30]. The S-nCCG algorithm is a decomposition method for two-stage robust optimization problem with continuous recourse [31]. In the master problem (MP), the middle maximization is replaced by critical uncertainty scenarios that are strategically identified by feasibility and optimality check sub-problems (SPs). Each SP requires an inner CCG loop.

To concisely present the proposed M-nCCG, we will use the general model of two-stage robust optimization problem with mixed-integer recourse, which is:

\{\begin{array}{l} \underset{γ \in Γ}{m i n} (p^{T} γ + \underset{ξ \in Ξ}{m a x} \underset{y, z}{m i n} (f (γ, ξ) + q^{T} y)) \\ s . t . A y + B z \leq b (γ, ξ) z i s b i n a r y \end{array}

(18)

where $γ$ is the vector of here-and-now variables; $ξ$ is the uncertainty vector; $y$ and $z$ are the wait-and-see variables; $A$ , $B$ , $b$ , $p$ , and $q$ are the coefficient matrices; and $f (γ, ξ)$ is the cost function. Note that (18) is an independent and pure math problem where the variables have no physical meanings.

1)　Outer CCG Loop

The MP in outer CCG loop (MP-outer) is:

\{\begin{array}{l} \underset{γ \in Γ, θ, (y_{i}, z_{i})}{m i n} (p^{T} γ + θ) \\ s . t . A y_{i} + B z_{i} \leq b (γ, ξ_{i}) \forall ξ_{i} \in Ξ^{†}, z_{i} i s b i n a r y \\ θ \geq f (γ, ξ_{i}) + q^{T} y_{i} \forall ξ_{i} \in Ξ^{†} \end{array}

(19)

where $Ξ^{†}$ contains the critical uncertainty scenarios that have been identified so far indexed by $i$ ; and $θ$ is an auxiliary variable.

Given $γ^{*}$ , there are two SPs. One is the feasibility check SP (SP1-outer), finding those uncertainty scenarios that cause the most severe infeasibility, i.e.,

\{\begin{array}{l} \underset{ξ \in Ξ}{m a x} \underset{μ \geq 0, y, z}{m i n} 1^{T} μ \\ s . t . A y + B z - μ \leq b (γ^{*}, ξ) z i s b i n a r y \end{array}

(20)

where $μ$ is a slack vector; and $1$ is an all-one column vector with compatible rows. If the optimal value of SP1-outer is zero, there is always a feasible $(y, z)$ for any $ξ \in Ξ$ ; otherwise, the optimal $ξ^{*}$ can cause infeasibility.

The other is the optimality check SP (SP2-outer), finding those uncertainty scenarios that cause the highest costs, i.e.,

\{\begin{array}{l} \underset{ξ \in Ξ}{m a x} \underset{y, z}{m i n} (f (γ, ξ) + q^{T} y) \\ s . t . A y + B z \leq b (γ, ξ) z i s b i n a r y \end{array}

(21)

2)　Inner CCG Loop for SP1-outer

SP1-outer and SP2-outer are intractable due to the binary $z$ , so the inner CCG loop is required. Considering SP1-outer first, we reformulate it as:

\{\begin{array}{l} \underset{ξ \in Ξ}{m a x} \underset{z i s b i n a r y}{m i n} \underset{μ \geq 0, y}{m i n} 1^{T} μ \\ s . t . A y - μ \leq b (γ^{*}, ξ) - B z : λ \end{array}

(22)

where $λ$ after the colon is the vector of dual variables. Dualizing the inner minimization leads to:

\{\begin{array}{l} \underset{ξ \in Ξ}{m a x} \underset{z i s b i n a r y}{m i n} \underset{λ}{m a x} λ^{T} (b (γ^{*}, ξ) - B z) \\ s . t . A^{T} λ = 0 - 1 \leq λ \leq 0 \end{array}

(23)

Note that (23) is a two-stage robust optimization problem with continuous recourse, which can be handled by a CCG loop whose MP is:

\{\begin{array}{l} \underset{ξ \in Ξ, θ, λ_{i}}{m a x} θ \\ \begin{array}{l} s . t . θ \leq λ_{i}^{T} (b (γ^{*}, ξ) - B z_{i}) \forall z_{i} \in 𝒵_{1}^{†} \\ A^{T} λ_{i} = 0 - 1 \leq λ_{i} \leq 0 \end{array} \end{array}

(24)

where $𝒵_{1}^{†}$ is the set of critical binary recourse actions that have been identified so far in the inner loop of SP1-outer. Given $ξ^{*}$ , the corresponding SP-inner1 is:

\{\begin{array}{l} \underset{z i s b i n a r y}{m i n} \underset{λ}{m a x} λ^{T} (b (γ^{*}, ξ^{*}) - B z) \\ s . t . A^{T} λ = 0 - 1 \leq λ \leq 0 \end{array}

(25)

In (22), the inner minimization is always feasible because there is a slack vector $μ$ . Hence, the feasible set of $λ$ in (23) is non-empty. Since the feasible set of $λ$ is independent from $ξ$ and $z$ , the inner maximization in (23) is also always feasible. Therefore, to solve (22) in the form of (23), we only need an SP-inner1. A feasibility check SP is unnecessary.

3)　Inner CCG Loop for SP2-outer

In the similar way, we can reformulate SP2-outer as:

\{\begin{array}{l} \underset{ξ \in Ξ}{m a x} \underset{z i s b i n a r y}{m i n} \underset{y}{m i n} (f (γ^{*}, ξ) + q^{T} y) \\ s . t . A y \leq b (γ^{*}, ξ) - B z : ζ \end{array}

(26)

where $ζ$ is the vector of dual variables. Using duality theory to minimize the inner, we have:

\{\begin{array}{l} \underset{ξ \in Ξ}{m a x} \underset{z i s b i n a r y}{m i n} \underset{ζ}{m a x} [ζ^{T} (b (γ^{*}, ξ) - B z) + f (γ^{*}, ξ)] \\ s . t . A^{T} ζ = q ζ \leq 0 \end{array}

(27)

Correspondingly, the MP-inner1 is:

\{\begin{array}{l} \underset{ξ \in Ξ, θ, ζ_{i}}{m a x} θ \\ s . t . θ \leq {ζ_{i}}^{T} (b (γ^{*}, ξ) - B z_{i}) + f (γ^{*}, ξ) \forall z_{i} \in 𝒵_{2}^{†} \\ A^{T} ζ_{i} = q ζ_{i} \leq 0 \end{array}

(28)

where $𝒵_{2}^{†}$ is the set of critical binary recourse actions that have been identified so far in the inner loop of SP2-outer.

The SP-inner2 is:

\{\begin{array}{l} \underset{z i s b i n a r y}{m i n} \underset{ζ}{m a x} [ζ^{T} (b (γ^{*}, ξ^{*}) - B z) + f (γ^{*}, ξ^{*})] \\ s . t . A^{T} ζ = q ζ \leq 0 \end{array}

(29)

For the same reason as in the last subsection, a feasibility check SP is not required here, as long as we can prove that the feasible set of dual vector $ζ$ is non-empty by finding at least one feasible primal vector $y$ . This can be achieved in the MP-outer by completing the feasibility check process before starting the optimality check process.

4)　M-nCCG Algorithm

The M-nCCG algorithm is finalized in Algorithm 1. Two remarks are given below.

Algorithm 1 : M-nCCG

Initialization: set $Ξ^{†} = \{ξ^{0}\}$ , $Z_{1}^{†} = \{z^{0}\}$ ,and $Z_{2}^{†} = \{z^{0}\}$ with arbitrary $ξ^{0}$ and z⁰; a tolerance $ε = 10^{- 6}$ ; a binary indicator $κ = 0$ ( $κ = 0 / 1$ means feasibility check is incomplete/complete).

Step 1: solve MP-outer and save the optimal $γ^{*}$ . The optimal value is $L B^{o u t e r}$ . If $κ = 1$ , go to Step 3.

Step 2: solve SP-outer by the following substeps.

Step 2-1: solve MP-inner1 and save the optimal solution $ξ^{*}$ and value $θ^{*}$ ; set the upper bound $U B^{i n n e r 1} = θ^{*}$ .

Step 2-2: solve SP-innerl and save the optimal solution $z^{*}$ . The optimal value is the lower bound, denoted as $L B^{i n n e r l}$ .

Step 2-3: if $U B^{i n n e r 1} - L B^{i n n e r 1} \geq ε$ , let $Z_{1}^{†} = Z_{1}^{†} ⋃ \{z^{*}\}$ and return to Step 2-1; if $U B^{i n n e r l} - L B^{i n n e r l} < ε$ and $L B^{i n n e r 1} = 0$ , set $κ = 1$ and go to Step 3; else if $U B^{i n n e r 1} - L B^{i n n e r l} < ε$ and $L B^{i n n e r l} > 0$ , let $Ξ^{†} = Ξ^{†} ⋃ \{ξ^{*}\}$ and return to Step 1.

Step 3: solve SP2-outer by the following substeps.

Step 3-1: solve MP-inner2 and save the optimal solution $ξ^{*}$ and value $θ^{*}$ ; set the upper bound $U B^{i n n e r 2} = θ^{*}$ .

Step 3-2: solve SP-inner2 and save the optimal solution $z^{*}$ . The optimal value is the lower bound, denoted as $L B^{i n n e r 2}$ .

Step 3-3: if $U B^{i n n e r 2} - L B^{i n n e r 2} \geq ε$ , let $Z_{2}^{†} = Z_{2}^{†} ⋃ \{z^{*}\}$ and return to Step 3-1; if $U B^{i n n e r 2} - L B^{i n n e r 2} < ε$ , set $U B^{o u t e r} = p^{T} γ^{*} + L B^{i n n e r 2}$ and go to Step 4.

Step 4: If $U B^{o u t e r} - L B^{o u t e r} < ε$ , terminate and report $γ^{*}$ ; else, $Ξ^{†} = Ξ^{†} ⋃ \{ξ^{*}\}$ and return to Step 1.

Output: $γ^{*}$ .

1) In the S-nCCG algorithm, the identified binary recourse actions ( $𝒵_{1}^{†}$ and $𝒵_{2}^{†}$ ) within each outer iteration are cleared when going to the next iteration. The proposed M-nCCG algorithm keeps them to give the next iteration a warm start, facilitating the convergence of inner loops.

A flowchart is given in Fig. 2, where the red boxes show the modification, which improves the performance by exploiting all the information in pervious iterations.

Fig. 2 Algorithm flowchart.

2) Both SP-inner1 and SP-inner2 are bi-level min-max problems. For each of them, with a fixed outer level, the inner level is a linear program. Therefore, dualizing the inner level will lead to an MILP, which can be efficiently solved by commercial solvers like CPLEX.

D. Implementation Procedure

Ahead of the day, problem (14) in the form of (15) is solved by Algorithm 1. The optimal solution ${\bar{φ}}^{*}$ includes UC strategy $φ^{*}$ and partially affine policy $(U_{t}^{*}, V_{t}^{*})$ , $\forall t \in T$ , and the former is applied and the latter is passed to real-time dispatch.

During period $t$ of real-time dispatch, we need to find $x_{t}^{*}$ , $y_{t}^{*}$ , and $z_{t}^{*}$ . By partially affine policy, we have $x_{t}^{*} = U_{t}^{*} ξ_{[t]} + V_{t}^{*}$ . Besides, we notice in (14) that all the optimization problems in the real-time dispatch are decoupled over periods. Therefore, we compute $y_{t}^{*}$ and $z_{t}^{*}$ by:

(y_{t}^{*}, z_{t}^{*}) = a r g \underset{y_{t}, z_{t}}{m i n} (r^{T} U_{t}^{*} ξ_{[t]} + s^{T} y_{t})

(30)

V. Case Studies

The proposed method is tested on the modified IEEE 118-bus system with 186 transmission lines. This system consists of 15 coal-fired units (2800 MW), 12 gas-fired units (1400 MW), and 6 wind farms (3200 MW). The system-wide load curve is drawn in Fig. 3(a), where the peak load is 4434.2 MW. The forecast sequence $\{ξ_{t}^{w, 0}\}$ , $t = 1 : T$ is drawn in Fig. 3(b), where the wind power is normalized with the installed capacity as base value. The uncertainty set is built with $δ^{w} = 0.3$ p.u.. We set $C_{l}^{S h e d} = 1000$ $/MWh and $C_{w}^{C u r t} = 15$ $/MWh.

Fig. 3 Load and wind data. (a) System-wide load. (b) Forecast sequence.

We choose 8 buses (substations) as candidates for load shedding. According to model (17), each of them has $N_{l}$ outlets in total, and $N_{l}^{*}$ outlets can be tripped. The load elsewhere must be satisfied. Table I presents the load shedding settings. The load at other buses must be satisfied. All the parameters and data used in this section are provided in an online archive [

32] for interested readers. Optimization problems are solved by CPLEX 12.6 on a laptop with Intel i5-8250U CPU and 8 GB memory.

TABLE I Load Shedding Settings

Index	Bus	Peak load (MW)	$N_{l}$	$N_{l}^{*}$
1	15	90	2	1
2	42	96	2	1
3	49	87	2	2
4	54	113	3	2
5	59	277	6	4
6	80	130	3	2
7	90	163	4	3
8	116	184	4	3

A. UC Results

By Algorithm 1, the MRUC strategy is obtained, including the on-off status of coal-fired and gas-fired units. The UC cost, i.e., the first term in (7), is $$ 1.942 \times 10^{5}$ .

All the coal-fired units are on over the entire horizon. One reason is that coal is cheaper than natural gas, so the coal-fired units keep running to serve the base load. The other reason is that coal-fired units are less flexible than gas-fired units and cannot switch the on-off status frequently.

Table II exhibits the UC of gas-fired units, where 1 means the unit is on and 0 means the unit is off.

TABLE II UC of Gas-fired Units

Period (hour)	UC
Period (hour)	1	2	3	4	5	6	7	8	9	10	11	12
1	0	0	0	0	0	0	0	1	1	0	0	0
2	0	0	0	1	0	0	0	1	1	0	0	0
3	0	0	0	0	0	0	0	1	1	0	0	0
4	0	0	0	0	0	0	0	1	1	0	0	0
5	0	0	0	0	1	0	0	1	1	0	0	0
6	0	1	0	0	1	0	0	1	1	0	0	0
7	0	1	0	0	1	0	0	1	1	1	0	0
8	0	0	0	0	1	0	0	1	1	1	0	0
9	0	1	0	0	1	0	0	1	1	0	0	0
10	1	1	1	1	1	1	1	1	1	1	1	0
11	1	1	1	1	1	1	1	1	1	1	1	0
12	1	1	1	1	1	1	1	1	1	1	1	1
13	1	1	1	1	1	1	1	1	1	1	1	1
14	1	1	1	0	1	1	1	1	1	1	1	0
15	1	1	1	0	1	1	0	1	1	1	1	0
16	1	1	1	0	1	1	0	1	1	1	1	0
17	1	1	1	1	1	1	1	1	1	1	1	0
18	1	1	0	1	1	1	0	1	1	1	1	0
19	1	1	1	1	1	1	1	1	1	1	1	1
20	1	1	1	1	1	1	1	1	1	1	1	1
21	1	1	1	1	1	1	1	1	1	1	1	1
22	1	1	1	1	1	1	1	1	1	1	1	0
23	0	1	0	0	1	0	0	1	1	0	0	0
24	0	1	0	0	1	0	0	1	1	1	0	0

1) Most gas-fired units are committed from hour 10 to hour 22, since the load within this time interval is relatively high. Meanwhile, units 8 and 9 keep running all day long to help serve the base load. Besides, unit 12 is deployed only at noon and in the evening, providing power support in peak hours.

2) The flexibility of gas-fired units can be observed from their quick status switches, e.g., unit 7 is turned on in hour 17, turned off soon in hour 18, and turned on again in hour 19. Such a capability enables to respond to the fast fluctuation of wind power.

These results show that gas-fired units are important for power supply, especially in terms of flexibility.

B. Dispatch Results in Worst-case Scenario

Besides UC strategy, Algorithm 1 also identifies the worst-case scenario, i.e., the one in $Ξ^{†}$ that causes the highest dispatch cost. The worst-case scenario here is exactly the lower bound in Fig. 3(b), where the wind energy is quite scarce.

The dispatch results with regard to this scenario are discussed below.

The coal fuel cost is $$ 9.602 \times 10^{5}$ in total, and the gas fuel cost is $$ 1.511 \times 10^{6}$ . There is no penalty of wind energy curtailment because wind energy is poor in the worst-case scenario. The total penalty of load shedding is $$ 5.581 \times 10^{5}$ , so the amount of load shedding is $558.1$ MWh.

Figure 4 shows the system-wide generation power curves of coal-fired units and gas-fired units. The output power of coal-fired units is stable and high. Since coal is cheaper than natural gas, the coal-fired units produce the majority of energy. Meanwhile, gas-fired units play a significant role in managing load peak and valley. The generation curve shares a similar trend with the load curve in Fig. 3(a). Around hour 4, the load is the lowest, so is the total output of gas-fired units. There exist two load peaks around hour 12 and hour 20 when the gas-fired units respond to enhance the output level.

Fig. 4 System-wide generation power curves of coal-fired units and gas-fired units.

Furthermore, Fig. 5 shows the load shedding over time. Load shedding happens in hours 19, 20, and 21 and the total amount is 558.1 MWh. The total load is $9.502 \times 10^{4}$ MWh over the entire horizon, so the load shedding ratio is about $0.6 %$ .

Fig. 5 Load shedding.

C. Sensitivity Analysis

The robustness of the proposed method is controlled by parameter $δ_{t}^{w}$ , which can be regarded as the maximum forecast error of wind power. Therefore, this subsection investigates the impact of $δ_{t}^{w}$ on UC and dispatch in the worst-case scenario.

The results of sensitivity analysis about $δ_{t}^{w}$ are gathered in Table III. With the increase of $δ_{t}^{w}$ , the wind energy in the worst-case scenario becomes less and less. Meanwhile, the possible fluctuation of wind power within the uncertainty set becomes more severe. Consequently, the UC cost goes up since the units are deployed for more time and with more status switches. The fuel consumption also increases to compensate the scarcity of wind energy, leading to a rising fuel cost. Especially, the gas fuel cost rises up by $$ 0.9150 \times 10^{6}$ when $δ_{t}^{w}$ changes from 0.1 to 0.5.

TABLE III Sensitivity Analysis About

δ_{t}^{w}

$δ_{t}^{w}$ (p.u.)	UC cost $(10^{5} $)$	Coal fuel cost $(10^{5} $)$	Gas fuel cost $(10^{6} $)$	Wind curtailment $(M W h)$	Load shedding $(M W h)$
0.1	1.358	9.214	0.880	141.90	0
0.2	1.699	9.504	1.233	84.73	204.2
0.3	1.942	9.602	1.511	0	558.1
0.4	2.044	9.627	1.707	0	595.8
0.5	2.101	9.635	1.795	0	595.8

Furthermore, the wind curtailment is zero when $δ_{t}^{w} \geq 0.3$ . A larger $δ_{t}^{w}$ implies that the wind energy in the worst-case scenario is rarer and thus can be fully utilized. On the contrary, load shedding is zero when $δ_{t}^{w} = 0.1$ , since the wind energy in the worst-case scenario is relatively rich.

D. Comparison with Fully Affine Policy

We propose the partially affine policy in (13) and the fully affine policy in (12) with $δ_{t}^{w} = 0.3$ . For both of them, binary recourse is optimized during each period.

We compute the sample average of real-time dispatch cost instead of the results regarding the worst-case scenario. To this end, we collect 200 per-unit wind samples (trajectories) from a real wind farm in Ningxia Province, China [

32]. Some samples fall within the uncertainty set

Ξ

and some not. For those within the uncertainty set, the real-time dispatch is always feasible. For those outside the uncertainty set, the infeasibility may occur.

Table IV summarizes the comparison results. The UC costs corresponding to these two policies are close.

TABLE IV Comparisons Between Partially and Fully Affine Policies

Item	UC cost $(10^{5} $)$	Average dispatch cost for 65 samples inside $Ξ$ $(10^{6} $)$	Number of feasible samples for 135 samples outside $Ξ$	Average dispatch cost for 135 samples outside $Ξ$ $(10^{6} $)$
Partially affine policy	1.666	1.593	86	1.541
Fully affine policy	1.682	1.827	69	1.800

In addition, we find that:

1) The proposed partially affine policy realizes a lower dispatch cost. According to (13), partially affine policy maintains the flexibility of non-state variables, which physically mean the outputs of gas-fired units; such flexibility can help find better real-time actions. With samples inside/outside the uncertainty, the dispatch cost of fully affine policy is significantly higher by $14.7 % / 16.8 %$ than that of partially affine policy. Therefore, the advantage of the proposed method over optimality is verified.

2) The flexibility maintained by partially affine policy favors the robustness. Among 135 samples outside the uncertainty, the real-time dispatch by partially affine policy is feasible with 86 samples, while the number becomes 69 by fully affine policy. The reason is that unleashing the flexibility of gas-fired units enlarges the feasible region of dispatch actions, so the resulting system is more robust to uncertain wind power. Therefore, the advantage of the proposed method over dispatch robustness is verified.

E. Comparison with S-nCCG: Computation Performance

According to Section IV-C, in each outer iteration, both the proposed M-nCCG algorithm and the S-nCCG algorithm solve the same MP (MP-outer). The difference is how they solve SP1-outer and SP2-outer by inner loops. After the $i$ ^th outer iteration, M-nCCG algorithm takes the identified binary recourse actions ( $𝒵_{1}^{†}$ and $𝒵_{2}^{†}$ ) to the $(i + 1)$ ^th outer iteration. The inner loops are accelerated by these actions as a warm start. S-nCCG renews $𝒵_{1}^{†}$ and $𝒵_{2}^{†}$ when stepping into the $(i + 1)$ ^th outer iteration.

In Fig. 6, M-nCCG and S-nCCG share the same outer convergence curve, which indicates that the algorithm terminates after 6 outer iterations. Taking the third outer iteration as an example, M-nCCG, which maintains the previously identified scenarios, takes only 2 iterations to terminate the inner loop of optimality check SP while S-nCCG takes 5 iterations. Therefore, M-nCCG can accelerate the convergence.

Fig. 6 Iteration of outer CCG loop and inner CCG loop. (a) Iteration of outer CCG loop. (b) Iteration of inner CCG loop.

The time to solve MP-outer is 103.6 s. To solve SP1-outer and SP2-outer, M-nCCG consumes 56.8 s while S-nCCG consumes 116.6 s. Therefore, M-nCCG reduces the total computation time by 27.2% compared with S-nCCG, which verifies the advantage over computation performance.

F. Comparison with State-of-the-art Methods

We compare the proposed method with two state-of-the-art methods mentioned in the literature review. The one in [

19], [20] denoted by M-1 uses the affine policy for continuous recourse and the piecewise constant policy for binary recourse. The piecewise constant policy relies on a pre-defined partitioning function. M-1 resorts to an MILP. The other one denoted by M-2 is reported in [21], which optimizes the partitioning function. M-2 employs a customized CCG algorithm. The total cost and computation time are compared, as shown in Table V. The total cost includes the UC cost and the sample average of dispatch cost.

TABLE V Comparisons with State-of-the-art Methods

Method	Total cost $(10^{6} $)$	Computation time (s)
Proposed	1.730	160.4
M-1	1.941	148.6
M-2	1.755	210.2

Regarding dispatch economy, M-1 achieves the highest total cost $$ 1.941 \times 10^{6}$ . The total cost of M-2 is 9.58% lower since it optimizes the partitioning function instead of using a pre-defined one like M-1. The total cost of the proposed method is the lowest. The reason is the proposed method does not impose any decision structure for non-state continuous and binary dispatch variables, whose flexibility is exploited and retained.

Regarding computation time, M-1 is the fastest and the proposed method has a close efficiency to M-1. M-2 consumes much more time than M-1 by 41.5%, since it entails solving a large-scale MP, which is established to improve the optimality. In summary, the proposed method balances well the optimality and efficiency.

VI. Conclusion

This paper investigates the UC problem considering discrete load shedding, which is formulated as a multi-stage adaptive robust optimization problem with mixed-integer recourse. Partially affine policy and two-stage reformulation address this problem in a tractable way. We conclude that the proposed partially affine policy outperforms fully affine policy in terms of dispatch economy and robustness, the M-nCCG algorithm can accelerate the convergence by warming up the inner loops, and the framework in this paper balances the optimality and efficiency compared with state-of-the-art methods.

The main weakness of the proposed method is that it can only handle linear problems. The exact power flow model is nonlinear and nonconvex. Besides, the power supply is affected by both balance and stability. The current formulation does not incorporate stability constraints. Future work will try to address these issues.

References

D. Bertsimas, E. Litvinov, X. Sun et al., “Adaptive robust optimization for the security constrained unit commitment problem,” IEEE Transactions on Power Systems, vol. 28, no. 1, pp. 52-63, Feb. 2013. [Baidu Scholar]

B. Hu and L. Wu, “Robust SCUC considering continuous/discrete uncertainties and quick-start units: a two-stage robust optimization with mixed-integer recourse,” IEEE Transactions on Power Systems, vol. 31, no. 2, pp. 1407-1419, Mar. 2016. [Baidu Scholar]

J. Zhang, Z. Chen, N. Zhang et al., “Frequency-constrained unit commitments with linear rules extracted from simulation results considering regulations from battery storage,” Journal of Modern Power Systems and Clean Energy, vol. 11, no. 4, pp. 1041-1052, Jul. 2023. [Baidu Scholar]

A. Lorca, X. Sun, E. Litvinov et al., “Multistage adaptive robust optimization for the unit commitment problem,” Operations Research, vol. 64, no. 1, pp. 32-51, Jan. 2016. [Baidu Scholar]

J. R. Birge and F. Louveaux, Introduction to Stochastic Programming, 2nd ed. New York, USA: Springer, 2011. [Baidu Scholar]

Z. Guo, W. Wei, L. Chen et al., “Distribution system operation with renewables and energy storage: a linear programming based multistage robust feasibility approach,” IEEE Transactions on Power Systems, vol. 37, no. 1, pp. 738-749, Jan. 2022. [Baidu Scholar]

A. Lorca and X. Sun, “Multistage robust unit commitment with dynamic uncertainty sets and energy storage,” IEEE Transactions on Power Systems, vol. 32, no. 3, pp. 1678-1688, May 2017. [Baidu Scholar]

Y. Zhou, M. Shahidehpour, Z. Wei et al., “Multistage robust lookahead unit commitment with probabilistic forecasting in multi-carrier energy systems,” IEEE Transactions on Sustainable Energy, vol. 12, no. 1, pp. 70-82, Jan. 2021. [Baidu Scholar]

N. G. Cobos, J. M. Arroyo, N. Alguacil et al., “Robust energy and reserve scheduling considering bulk energy storage units and wind uncertainty,” IEEE Transactions on Power Systems, vol. 33, no. 5, pp. 5206-5216, Sept. 2018. [Baidu Scholar]

Y. Zhou, Q. Zhai, and L. Wu, “Multistage transmission-constrained unit commitment with renewable energy and energy storage: implicit and explicit decision methods,” IEEE Transactions Sustainable Energy, vol. 12, no. 2, pp. 1032-1043, Apr. 2021. [Baidu Scholar]

X. Zheng, M. E. Khodayar, J. Wang et al., “Distributionally robust multistage dispatch with discrete recourse of energy storage systems,” IEEE Transactions on Power Systems, doi: 10.1109/TPWRS.2024.3369664 [Baidu Scholar]

A. Georghiou, A. Tsoukalas, and W. Wiesemann, “Robust dual dynamic programming,” Operations Research, vol. 67, no. 3, pp. 813-830, May 2019. [Baidu Scholar]

Y. Shi, S. Dong, C. Guo et al., “Enhancing the flexibility of storage integrated power system by multi-stage robust dispatch,” IEEE Transactions on Power Systems, vol. 36, no. 3, pp. 2314-2322, May 2021. [Baidu Scholar]

E. Nasrolahpour, J. Kazempour, H. Zareipour et al., “A bilevel model for participation of a storage system in energy and reserve markets,” IEEE Transactions on Sustainable Energy, vol. 9, no. 2, pp. 582-598, Apr. 2018. [Baidu Scholar]

X. Xu, H. Zhang, C. Li et al., “Optimization of the event-driven emergency load-shedding considering transient security and stability constraints,” IEEE Transactions on Power Systems, vol. 32, no. 4, pp. 2581-2592, Jul. 2017. [Baidu Scholar]

H. Qiu, W. Gu, W. Sheng et al., “Resilience-oriented multistage scheduling for power grids considering nonanticipativity under tropical cyclones,” IEEE Transactions on Power Systems, vol. 38, no. 4, pp. 3254-3267, Jul. 2023. [Baidu Scholar]

D. Bertsimas and C. Caramanis, “Finite adaptability in multistage linear optimization,” IEEE Transactions on Automatic Control, vol. 55, no. 12, pp. 2751-2767, Dec. 2010. [Baidu Scholar]

A. Subramanyam, C. E. Gounaris, and W. Wiesemann, “K-adaptability in two-stage mixed-integer robust optimization,” Mathematical Programming Computation, vol. 12, pp. 193-224, Dec. 2020. [Baidu Scholar]

D. Bertsimas and A. Georghiou, “Binary decision rules for multistage adaptive mixed-integer optimization,” Mathematical Programming, vol. 167, pp. 395-433, Jan. 2018. [Baidu Scholar]

H. Qiu and H. B. Gooi, “A unified MILP solution framework for adaptive robust scheduling problems with mixed-integer recourse objective,” IEEE Transactions on Power Systems, vol. 38, no. 1, pp. 952-955, Jan. 2023. [Baidu Scholar]

Z. Zhong, N. Fan, and L. Wu, “Multistage robust optimization for the day-ahead scheduling of hybrid thermal-hydro-wind-solar systems,” Journal of Global Optimization, vol. 27, pp. 1-36, Nov. 2023. [Baidu Scholar]

Q. Chen, P. Zou, C. Wu et al., “A Nash-Cournot approach to assessing flexible ramping products,” Applied Energy, vol. 206, pp. 42-50, May 2017. [Baidu Scholar]

C. Wang, R. Gao, W. Wei et al., “Risk-based distributionally robust optimal gas-power flow with Wasserstein distance,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 2190-2204, May 2019. [Baidu Scholar]

A. Belderbos, T. Valkaert, K. Bruninx et al., “Facilitating renewables and power-to-gas via integrated electrical power-gas system scheduling,” Applied Energy, vol. 275, p. 115082, Jan. 2020. [Baidu Scholar]

P. Xiong and C. Singh, “A distributional interpretation of uncertainty sets in unit commitment under uncertain wind power,” IEEE Transactions on Sustainable Energy, vol. 10, no. 1, pp. 149-157, Jan. 2019. [Baidu Scholar]

C. Zhao and Y. Guan, “Unified stochastic and robust unit commitment,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 3353-3361, Aug. 2013. [Baidu Scholar]

A. Ben-Tal, A. Goryashko, E. Guslitzer et al., “Adjustable robust solutions of uncertain linear programs,” Mathematical Programming, vol. 99, pp. 351-376, Jun. 2004. [Baidu Scholar]

D. Bertsimas and A. Georghiou, “Binary decision rules for multistage adaptive mixed-integer optimization,” Mathematical Programming, vol. 167, pp. 395-433, Jul. 2018. [Baidu Scholar]

C. Ning and F. You, “A transformation-proximal bundle algorithm for multistage adaptive robust optimization and application to constrained robust optimal control,” Automatica, vol. 113, p. 108802, Dec. 2020. [Baidu Scholar]

L. Zhao and B. Zeng. (2012, Jan.). An exact algorithm for two-stage robust optimization with mixed integer recourse problems. [Online]. Available: https://optimization-online.org/wpcontent/uploads/2012/01/3310.pdf. [Baidu Scholar]

B. Zeng and L. Zhao, “Solving two-stage robust optimization problems using a column-and-constraint generation method,” Operations Research Letters, vol. 41, pp. 457-461, Sept. 2013. [Baidu Scholar]

Z. Guo. (2023, Feb.). Data MRUC integer recourse. [Online]. Available: https://github.com/ZhongjieGuo/Papers. [Baidu Scholar]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher