A Continuous Operating Envelope for Managing Intra-interval Fluctuations: Modeling and Solution

Menghan Zhang; Zhifang Yang; Juan Yu; Wenyuan Li

网刊加载中。。。

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

确定继续浏览么?

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

A Continuous Operating Envelope for Managing Intra-interval Fluctuations: Modeling and Solution PDF

- ORCID：
Menghan Zhang
✉
- ORCID：
Zhifang Yang
✉
- ORCID：
Juan Yu
✉
- ORCID：
Wenyuan Li
✉

State Key Laboratory of Power Transmission Equipment Technology, School of Electrical Engineering, Chongqing University, Chongqing 400044, China

Updated：2025-03-26

DOI：10.35833/MPCE.2024.000636

OUTLINE

Abstract

Maintaining a continuous power balance is crucial for ensuring operational feasibility in power systems. However, due to forecasting difficulties and computational limitations, economic dispatch often relies on discrete interval horizons, which fail to guarantee feasibility within each interval. This paper introduces the concept of a continuous operating envelope for managing intra-interval fluctuations, delineating the range within which fluctuations remain manageable. We propose a parametric programming model to construct the envelope, represented as a polytope that accounts for both timescale and fluctuation dimensions. To address the computational challenges inherent in the parametric programming model, we develop a fast solution method to provide an approximated polytope. The approximated polytope, initially derived from lower-dimensional projections, represents a subset of the exact polytope that ensures operational feasibility. Additionally, we apply a polytope expansion strategy in the original dimensions to refine the approximated polytope, bringing the approximation closer to the exact polytope. Case studies on an illustrative 5-bus and a utility-scale 661-bus system demonstrate that the method effectively and stably provides a continuous operating envelope, particularly for high-dimensional problems.

Keywords

Continuous operating envelope; economic dispatch; intra-interval fluctuation; feasibility; uncertainty; polytope; parametric programming; ancillary service

I. Introduction

THE primary duty of a power system is to guarantee just-in-time production and transmission of bulk power to meet just-in-time customer demand [

1]. Although maintaining power balance over a continuous time horizon would be ideal, it is often impractical due to forecasting difficulties and computational limitations. Instead, system operators commonly employ a discrete interval horizon [2]. For example, in China, a 15-min horizon is used for the day-ahead dispatch. In the US, the independent system operators (ISOs) use an hourly horizon for the day-ahead dispatch and a 5-min horizon for the real-time operations.

To mitigate the limitations of discrete interval horizons, ISOs can rely on ancillary services to manage uncertainties within these intervals. However, the operational feasibility of dispatch under intra-interval fluctuations is not guaranteed, which may lead to reliability issues during real-time operations [

3]. With the increasing integration of renewable energy sources into the power system, their inherent intermittency and volatility heighten the need for timely assessments of operational feasibility under intra-interval fluctuations. Establishing a clear operating range, termed as the continuous operating envelope in this paper, is crucial for evaluating whether upcoming fluctuations threaten system operations. The continuous operating envelope provides clearer physical insights by explicitly defining the allowable range of intra-interval fluctuations, ensuring operational feasibility.

In practice, ancillary services are widely used to mitigate the impact of intra-interval fluctuations. For example, midcontinent ISO (MISO) and California ISO (CAISO) have developed flexible ramping products [

4], while Electric Reliability Council of Texas (ERCOT) has offered a fast regulation response [5]. Additionally, CAISO has introduced the concept of regulation mileage as an ancillary service [6]. These services improve the system ability to respond to variable demand and generation. However, these services lack a well-defined range of fluctuations that can be effectively managed.

Notably, fluctuations in the interval generally have a tighter allowable range compared with those occurring later. This phenomenon reflects temporal variance in the system ability to manage fluctuations. Relying solely on discrete interval horizons fails to capture this variance adequately, resulting in reactive rather than proactive fluctuation management. For example, the transitioning from day-ahead hourly results to real-time 5-min operations often requires reactive measures, such as curtailments or additional adjustments, in response to violations involving ramping or transmission issues [

7].

The high-resolution formulation offers an alternative method for managing intra-interval fluctuations. In this formulation, the number of discrete intervals throughout the decision period is increased to provide a finer granularity of dispatch and to better capture fluctuations [

8]. Intervals are set to be much denser, allowing them to approximate a continuous curve [9]. This method improves the power balance by considering all expected values within the forecasted range. However, it fails to establish an allowable range of fluctuations. Moreover, this method significantly increases computational burden at higher resolutions and the feasibility in shorter dispatch intervals remains uncertain [10]. Therefore, the high-resolution formulation is less commonly used for the day-ahead dispatch decisions.

Stochastic optimization offers another method for managing fluctuations. It models uncertainties through specific scenarios [

11] or uncertainty sets [12], which are determined by pre-specified probability distributions [13]. These distributions theoretically define the allowable range of fluctuations. However, these methods focus exclusively on power balance uncertainties at specific time steps within discrete interval horizons. They fail to provide a comprehensive description for managing intra-interval fluctuations, especially given that the system ability to manage fluctuations varies over time.

Moreover, region-based methods have been developed to delineate the allowable range of fluctuations caused by variable demand and generation [

14], [15]. The security region concept is first introduced for transmission networks [16], [17] and is later extended to distribution networks [18], [19]. The methods of active-power steady-state security region provide insights into ramping constraints and describe the corresponding polytope for each time step [20]. Recent research works have further extended these methods to incorporate various stability constraints under different operational conditions [21], [22]. However, security regions cannot explicitly represent polytope hyperplanes over time, thus requiring extensive sampling to approximate continuous outcomes. Additionally, security regions are determined solely by network topology, neglecting operational conditions resulting from dispatch decisions [23], [24].

In market-based power systems, day-ahead dispatch results require additional guidances to manage operational risks from intra-interval fluctuations and to adjust dispatch outcomes for active power flow. For instance, out-of-market corrections in the US market [

25] and short-term adequacy assessment in the EU market [26] are used to adjust generation and network constraints based on the updated operational conditions for active power flow.

In summary, managing intra-interval fluctuations remains a significant challenge under a discrete interval horizon. Current research methods lack a well-defined allowable range to ensure feasibility within each interval, which may neglect operational risks. Therefore, a new method is needed to define an allowable range for intra-interval fluctuations while considering temporal variance in the system ability over time.

In this paper, we introduce the concept of a continuous operating envelope, which supports fluctuation management by accounting for temporal variations in the system ability (represented by ramping capacities in this paper). The envelope is represented by a theoretical polytope, which delineates the allowable range of fluctuations that the system can sustain without compromising reliability. The contributions are twofold.

1) We present a parametric programming model for theoretically constructing the continuous operating envelope, capturing the system ability to manage intra-interval fluctuations. The model uses a parametric timescale to describe time-varying constraints and treats fluctuations as independent parameters, resulting in a polytope that incorporates both timescale and fluctuations in the parameter space. Furthermore, we provide a theoretical analysis of the computational complexity of the model.

2) We develop a fast solution method to stably construct a continuous operating envelope. To reduce computational complexity, this method explores multiple lower-dimensional projections and then algebraically reconstructs the approximated polytope in the original dimensions. As a subset of the exact solution, the approximated polytope ensures full feasibility within its interior, making it suitable to be used as the continuous operating envelope. Furthermore, we propose a polytope expansion strategy to refine the approximated polytope, bringing it closer to the exact polytope.

The continuous operating envelope serves as an analytical tool to determine whether fluctuations are feasible and whether they pose a threat to system operations. Case studies include a illustrative 5-bus system to explain the basic concept and a utility-scale 661-bus system to verify the scalability. Results show the proposed method accurately constructs the envelope for the illustrative 5-bus system and successfully generates the approximated polytope as the envelope for the utility-scale 661-bus system. The results in the utility-scale 661-bus system also validate the effectiveness of the proposed method for high-dimensional problems.

II. Mathematical Formulation of Continuous Operating Envelope

A. Feasibility Examination for Intra-interval Fluctuations

Feasibility examination for intra-interval fluctuations plays a crucial role in proactively identifying and mitigating operational risks. Day-ahead dispatch provides hourly results, and feasibility examinations can be conducted at any time step within these intervals. It helps ensure the feasibility in response to forecasted intra-interval demand and renewable energy output fluctuations. If operations are found to be infeasible, operators can promptly mitigate risks through curtailments or other corrective measures.

Assume the feasibility examination is conducted at a specific time step t. The time step t lies within the interval $𝒯$ , defined by $t^{S} \leq t \leq t^{E}$ , where $t^{S}$ and $t^{E}$ are the start and end of the interval, respectively. To conduct the feasibility examination, the direct current (DC) optimal power flow (OPF) model is employed, which includes:

1)　Objective function

\underset{G (t)}{m i n} z = H_{1} G (t) + H_{2} e_{G}^{T}

(1)

2)　Power balance constraint

e_{G}^{} G (t) + e_{R} P_{R}^{} (t) = e_{D}^{} P_{D}^{} (t)

(2)

3)　Transmission constraints

P_{L}^{m i n} \leq P_{L}^{} (t) \leq P_{L}^{m a x}

(3)

P_{L}^{} (t) = Γ (Ζ_{G}^{} G (t) + Ζ_{R}^{} P_{R}^{} (t) - Ζ_{D}^{} P_{D}^{} (t))

(4)

4)　Ramping constraint

(t - t^{S}) R_{D} \leq G (t) - G (t^{S}) \leq (t - t^{S}) R_{U}

(5)

5)　Generation constraint

G_{}^{m i n} \leq G (t) \leq G_{}^{m a x}

(6)

where $H_{1}$ and $H_{2}$ are the row vectors representing the marginal production costs and no-load costs of online units, respectively; $z$ is the objective function value; $G (t)$ is the column vector of decision variables representing the unit output (dispatchable resources) at time step t; $P_{R} (t)$ and $P_{D} (t)$ are the column vectors of the forecasted values representing the fluctuations in intra-interval renewable energy output and demand at time step $t$ , respectively; $e_{G}$ , $e_{R}$ , and $e_{D}$ are the row vectors filled with “all ones” associated with $G (t)$ , $P_{R} (t)$ , and $P_{D} (t)$ , respectively; $P_{L} (t)$ is the power flow at time step $t$ ; $Z_{G}$ , $Z_{R}$ , and $Z_{D}$ are the corresponding incident matrices; $Γ$ is the matrix of power transfer distribution factor; $R_{D}$ and $R_{U}$ are the ramping capacities per-unit time of dispatchable resources; superscripts max and min denote the maximum and minimum limits of variables and parameters, respectively; the superscript T is the transpose of a matrix; and $G (t^{S})$ is the unit output at time $t^{S}$ based on the given discrete dispatch decisions, such as the day-ahead hourly results.

In our analysis, temporal variations in ramping capacities at time step $t$ play a crucial role in determining the system ability to manage fluctuations, as these capacities typically vary across time steps [

27]. To emphasize the impact of temporal constraints, we express

R_{D}

and

R_{U}

as per-unit time values, which makes time step t a significant factor in (5).

$P_{R} (t)$ and $P_{D} (t)$ are non-dispatchable resources. They can be represented as forecasted functions of time step $t$ . At any node, $P_{R} (t)$ and $P_{D} (t)$ can be treated as fluctuations, while other nodes are considered non-fluctuating nodes. The model in (1)-(6) conducts a feasibility examination to ensure operational feasibility for these forecasted fluctuations. To facilitate the analysis, $P_{R} (t)$ and $P_{D} (t)$ are combined into a single vector $F (t)$ . This vector represents deviations at specific nodes resulting from the forecasted fluctuations at time step $t$ and is expressed as:

F (t) = Z_{R} (P_{R} (t) - P_{R} (t^{S})) + Z_{D} (P_{D} (t) - P_{D} (t^{S}))

(7)

where $P_{R} (t^{S})$ and $P_{D} (t^{S})$ are the renewable energy output and demand at the start of the interval $(t = t^{S})$ based on the given discrete dispatch decisions, respectively. If $P_{R} (t)$ and $P_{D} (t)$ deviate from $P_{R} (t^{S})$ and $P_{D} (t^{S})$ , respectively, this indicates that fluctuations are present at fluctuating nodes. $P_{R} (t)$ and $P_{D} (t)$ can fluctuate throughout the interval $𝒯$ . In (7), $F (t)$ represents all nodes, but for clarity, the following analysis will focus solely on the fluctuating nodes, excluding the non-fluctuating ones. To simplify visualization, we display the fluctuations $F (t)$ as $P_{R} (t)$ and $P_{D} (t)$ in the subsequent figures, without subtracting the constant values $P_{R} (t^{S})$ and $P_{D} (t^{S})$ at $t = t^{S}$ .

To focus on the key aspects, the feasibility examination model in (1)-(6) can be reformulated as the following general linear optimization problem in (8) with $G (t)$ , $F (t)$ , and the given time step $t$ . Let n denote the dimension of the decision variables, and $G (t)$ is an $n \times 1$ decision variable vector. Let $m$ denote the dimension of the fluctuation parameters, and $F (t)$ is an $m \times 1$ parameter vector representing $m$ fluctuating nodes. Let $p$ denote the dimension of the constraints in (8). A is a $p \times n$ constant matrix. $B$ is a $p \times 1$ column vector. $C$ is a $p \times m$ constant matrix. $D$ is a $p \times 1$ column vector. $𝒢$ is the polytope of $G (t)$ . We can obtain:

\{\begin{array}{l} \underset{G (t)}{m i n} z = H_{1} G (t) + H_{2} e_{G}^{T} \\ s . t . A G (t) \leq B + C F (t) + t D G (t) \in 𝒢 \subseteq ℝ^{n} \end{array}

(8)

In this optimization problem, $F (t)$ is treated as a constant vector with the forecasted values. By collecting all feasible $F (t)$ samples in (8), we define a closed range that ensures feasibility for each time step $t$ . Figure 1 visualizes these findings using 2-dimensional $F (t)$ samples ( $F_{1} (t)$ and $F_{2} (t)$ ) and the time step $t$ . In Fig. 1, green points represent infeasible $F (t)$ samples, while the blue points represent feasible ones. As more $F (t)$ samples are found to be feasible, the points gradually approximate a closed blue line, indicating the allowable range of fluctuations at time step $t$ . However, achieving the theoretical blue line would require an infinite number of $F (t)$ samples.

Fig. 1 Feasibility examination for intra-interval fluctuations. (a) Closed allowable range for a single time step. (b) Allowable range for entire interval.

Moving beyond a single time step, we extend the analysis by continuously taking samples of $F (t)$ across the time steps within the interval $𝒯$ . The red point represents the initial $F (t^{S})$ at the start of the interval $(t = t^{S})$ based on the given discrete dispatch decision. This process gradually accumulates a series of blue points, eventually shaping a polytope. With an infinite number of samples, this polytope can be precisely defined. As shown in Fig. 1(b), the boundaries of the polytope, marked in blue, delineate the overall allowable range of intra-interval fluctuations.

In the rest of this paper, this polytope, which delineates the allowable range of intra-interval fluctuations, is defined as the continuous operating envelope. However, the model in (8) cannot analytically provide the continuous operating envelope from a finite number of samples. While the DC OPF method is effective for discrete analysis, it requires a more efficient and theoretical method for the continuous examination. The new method should minimize the repeated calculations, eliminate the need for infinite samples, and provide a more accurate and efficient feasibility examination.

B. Parametric Programming for Intra-interval Fluctuations

The continuous operating envelope physically represents the system ability to manage intra-interval fluctuations and defines the extreme operational boundaries. Revisiting the optimization problem in (8), the feasibility and the value of the objective function depend on given values of $F (t)$ and $t$ .

By treating $F (t)$ and $t$ as varying parameters, we formulate a parametric programming model for the continuous operating envelope. In this model, all fluctuations are considered independent parameters and are no longer linked to $t$ . To distinguish them from the original optimization problem, the parametric fluctuation vector of dimension $m \times 1$ is denoted by $F^{P}$ , and the superscript P denotes the parametric programming.

\{\begin{array}{l} \underset{G (F^{P}, t)}{m i n} z = H_{1} G (F^{P}, t) + H_{2} e_{G}^{T} \\ s . t . A G (F^{P}, t) \leq B + C F^{P} + t D \\ G (F^{P}, t) \in 𝒢 \subseteq ℝ^{n}, F^{P} \in ℱ \subseteq ℝ^{m}, t \in 𝒯 \subseteq ℝ^{1} \end{array}

(9)

where $ℱ$ is the polytope of $F^{P}$ .

We combine $F^{P}$ and $t$ into a single vector $θ$ of dimension $(m + 1) \times 1$ . Thus, $G (F^{P}, t)$ can be written as $G (θ)$ .

θ : = [\begin{array}{l} F^{P} \\ t \end{array}] θ \in Θ \subseteq ℝ^{m + 1}

(10)

where $Θ$ is the polytope of $θ$ .

In (9), $F^{P}$ at the fluctuating nodes is treated as a variable, meaning its specific values do not need to be known in advance, unlike that in (8). For clarity, in the following text, $G (t)$ and $F (t)$ are used in the optimization problem in (8), while $G (θ)$ and $F^{P}$ are used in the parametric programming model.

In parametric solutions, $G (θ)$ is characterized by $θ$ . Once $θ$ is fully specified, the model in (9) becomes equivalent to the optimization problem in (8). For example, the samples of $F (t)$ and t, which are known when executing the feasibility examination, can be regarded as a specific set-point $θ^{*}$ of $θ$ . The model in (9) provides the decision variable vector $G (θ^{*})$ , which is identical to the optimal result $G (t)$ from the model in (8) with $F (t)$ and $t$ . Therefore, all feasible set-points in (8) belong to $Θ$ , and $Θ$ serves as the theoretical polytope representation of the continuous operating envelope.

The process is outlined for determining $Θ$ using the model in (9). By selecting a $θ^{*}$ , $G (θ^{*})$ , Lagrange multiplier $λ (θ^{*})$ , and objective value $z (θ^{*})$ can be obtained from the model in (9). At this point, the model in (9) is reduced to an optimization problem. According to sensitivity analysis theory [

28], after optimization, the values in the neighborhood of

θ^{*}

can be expressed as affine functions of the varying

θ

\{\begin{array}{l} G^{P} (θ) = h_{G} (θ) \\ z^{P} (θ) = h_{z} (θ) \\ λ^{P} (θ) = h_{λ} (θ) \end{array}

(11)

where $h_{G}$ , $h_{z}$ , and $h_{λ}$ are the affine functions in parametric solutions of $G (θ)$ , $z (θ)$ , and $λ (θ)$ , respectively.

By introducing the affine functions $G^{P} (θ)$ into the constraint in (9), the constraints are classified as active and inactive constraints, which become parametric constraints.

A G^{P} (θ) \leq B + C F^{P} + t D \Rightarrow \{\begin{matrix} h_{J} (F^{P}, t) = h_{J} (θ) = 0 \\ h_{I} (F^{P}, t) = h_{I} (θ) \leq 0 \end{matrix}

(12)

where the subscripts $J$ and $I$ are the active and inactive constraint sets, respectively; and $h_{J}$ and $h_{I}$ are the affine functions for active and inactive constraint sets, respectively.

$G^{P} (θ)$ bounded by the same inactive constraints can be gathered in one region $Ω$ as:

Ω : = {θ \in ℝ^{m + 1} : h_{I} (θ) \leq 0}

(13)

where the region $Ω$ is a polyhedron that represents the inactive constraint set, and it is defined by the parametric constraints.

The feasibility condition is ensured by substituting the $G^{P} (θ)$ into the inactive constraints given by $Ω$ . Additionally, the optimality condition is given by $λ^{P} (θ) \geq 0$ , which adheres to the Karush-Kuhn-Tucker optimality condition. In accordance with parametric programming theory [

29], [30] and taking into account the parameter bounds, the critical region

ℋ

is defined by:

ℋ : = {θ \in ℝ^{m + 1} : h_{I} (θ) \leq 0, λ^{P} (θ) \geq 0, θ \leq θ^{m a x}}

(14)

In the above analysis, $h_{G}$ , $h_{λ}$ , $h_{z}$ , $h_{J}$ , and $h_{I}$ are calculated at $θ^{*}$ . The corresponding $ℋ$ refers to a region in the parameter space, where the structure of the affine functions remains unchanged.

In order to distinguish each critical region and its associated affine functions, we introduce the subscript $i \in ℕ^{+}$ as the region index, involving $ℋ_{i}$ , $G_{i}^{P} (θ) = h_{G, i} (θ)$ , $λ_{i}^{P} (θ) = h_{λ, i} (θ)$ , $z_{i}^{P} = h_{z, i} (θ)$ , $h_{J, i} (θ)$ , and $h_{I, i} (θ)$ . In this paper, a collection of objects, e.g., ${ℋ_{1}, \dots, ℋ_{i}, \dots, ℋ_{ℐ}}$ , is denoted by ${ℋ_{i}}^{ℐ}$ , where the superscript $ℐ$ is the total number of critical regions.

Once all the critical regions are identified, ${ℋ_{i}}^{ℐ}$ is collected and termed a polyhedral partition of $Θ$ . $Θ$ , which is defined by the outermost boundaries of ${ℋ_{i}}^{ℐ}$ , represents the continuous operating envelope.

Θ : = b d ({ℋ_{i}}^{ℐ})

(15)

where $b d$ is the outermost facet of a polyhedral partition.

As shown in (13)-(15), $Θ$ is expressed in terms of its halfspace representation (H-rep), which is the intersection of a finite number of halfspaces. An equivalent representation is the vertex representation (V-rep). The polytope is defined as $Θ : = c o n v ({𝒱_{i}}^{ℐ})$ , where ${𝒱_{i}}^{ℐ} = {𝒱_{1}, \dots, 𝒱_{i}, \dots, 𝒱_{ℐ}}$ , and $𝒱_{i}$ is the set of vertices corresponding to the critical region $ℋ_{i}$ .

To illustrate the process for determining $Θ$ , we visualize the 2-dimensional $F^{P}$ ( $F_{1}^{P}$ and $F_{2}^{P}$ ) and $t$ . In Fig. 2(a), we initially identify a critical region $ℋ_{i}$ (green region) at a specific $θ^{*}$ (blue point) from (9), which executes a single optimization problem with a given $θ^{*}$ . The parametric programming model then explores beyond this critical region to identify new critical regions at different set-points. This exploration continues until no additional critical regions emerge. Figure 2(b) displays the outcome: five critical regions ( $ℐ = 5$ ) identified by using at least five feasible set-points, each marked by a unique color. Finally, $Θ$ is defined by merging the five regions ${ℋ_{i}}^{5}$ according to (15).

Fig. 2 Continuous operating envelope from parametric programming model. (a) A critical region $ℋ_{i}$ with $θ^{*}$ . (b) $Θ$ in parameter space.

C. Computational Complexity of Parametric Programming

The advantage of the parametric programming model is its ability to represent decision variables as affine functions of parameters, eliminating the need for exhaustive sampling in (8). For the envelope problem, the model in (9) offers an efficient alternative to exhaustive sampling of $F (t)$ and $t$ . Using $F^{P}$ and $t$ , it generates $G_{i}^{P} (θ)$ for $ℋ_{i}$ and constructs ${ℋ_{i}}^{ℐ}$ . Thorough exploration of all critical regions is crucial for accurately defining $Θ$ .

The computational complexity of the model in (9) is influenced by both the number of active constraints and the dimensionality of the parameter space. The computational complexity is primarily determined by the number of critical regions. Assume there is a given set of $p$ constraints, and at any set-point in the parameter space, a maximum of $q$ constraints can be active. Thus, the number of possible combinations of active constraints, denoted by $η$ , is equal to (16) in the worst case.

\{η = \sum_{l = 0}^{q} (\begin{matrix} p \\ l \end{matrix}) (\begin{matrix} p \\ l \end{matrix}) = \frac{p!}{(p - l)! l!}

(16)

According to [

31], the number of critical regions

α

is bounded in the worst case.

α \leq \sum_{g = 0}^{η - 1} g! p^{g}

(17)

where g can be regarded as the search tree level.

As the number of parameters increases, the dimensionality of the parameter space expands, which further increases the computational complexity. For a closed polyhedron of $ℋ_{i}$ in $(m + 1)$ -dimensional parameter space, at least $m + 2$ halfspaces are required. According to [

32], for a polyhedron with

m + 2

halfspaces, the size of the unexplored regions

β

is calculated as:

β = \sum_{l = 0}^{m + 1} (\begin{matrix} m + 2 \\ l \end{matrix}) - 1

(18)

Typically, a closed polyhedron in an $(m + 1)$ -dimensional parameter space has $m + ε$ halfspaces, where $ε \geq 2$ . As $ε$ increases, the size of the unexplored regions increases. For simplicity, we assume that $ε = 2$ for all polyhedra. The required number of set-points to identify all critical regions is:

μ \geq α β

(19)

To identify the entire ${ℋ_{i}}^{ℐ}$ , at least $α$ feasible set-points and $μ$ set-points are required. The total number $μ$ depends on the method used to select set-points but is subject to a lower bound. In high-dimensional parameter spaces, this lower bound grows exponentially. Therefore, addressing the inherent computational challenges is crucial to improving the applicability of the parametric programming model. For the envelope problem, this paper proposes a fast solution method to calculate the polytope, which will be discussed in the next section.

III. Fast Solution Method of Continuous Operating Envelope

A. Main Idea of Fast Solution Method

The primary aim of this paper is to develop a theoretical method for assessing the feasibility of intra-interval fluctuations. In Section II, we introduce the concept of a continuous operating envelope, represented by $Θ$ . According to (15), $Θ$ is defined by the outermost boundaries of ${ℋ_{i}}^{ℐ}$ , and constructing $Θ$ does not require identifying all critical regions. Instead, it can be constructed by focusing on a smaller subset of critical regions that encompass the outermost ones, which significantly reduces computational complexity by excluding certain critical regions.

One major challenge is discerning which critical regions are essential and which are less important. To address this, we analyze the characteristics of fluctuations. In high-dimensional parameter spaces, identifying critical regions, where multiple fluctuations interact simultaneously, requires a significant number of set-points. These critical regions, however, are absent in lower-dimensional parameter spaces.

The fast solution method prioritizes exploring critical regions in low-dimensional projections. Since ramping capacities are defined by the timescale, we focus on 2-dimensional projections that include both a single fluctuation and the timescale. Through algebraic manipulation, critical regions in these 2-dimensional projections are transformed back into their original dimensions to form an approximated polytope. Additionally, a polytope expansion strategy is introduced to refine the approximated polytope by subdividing unexplored regions.

Mathematically, the parametric programming model in high-dimensional spaces often faces the curse of dimensionality. By analyzing high-dimensional results through low-dimensional projections, we offer a stable method for obtaining the approximated polytope, which reduces the initial complexity and lays the foundation for further refinement. The polytope expansion strategy improves the approximated polytope, aiming for greater accuracy and possibly achieving the exact polytope. This two-step method provides a stable envelope for high-dimensional problems, with the approximated polytope followed by more accurate refinement.

In summary, the proposed fast solution method stably constructs an approximated polytope, serving as a continuous operating envelope. While ensuring operational feasibility, the approximated polytope introduces a degree of conservatism since it is a subset of the accurate polytope. The polytope expansion strategy balances computational efficiency with accuracy, which aims to progressively refine the polytope.

B. Low-dimensional Projections for Single Fluctuation

As the dimensionality of the parameter space increases, assigning suitable set-points to identify critical regions becomes increasingly difficult. In an $(m + 1)$ -dimensional space, each $θ^{*}$ for $θ$ corresponds to an $(m + 1)$ -dimensional point.

θ : = [\begin{array}{l} F^{P} \\ t \end{array}] = [F_{1}^{P}, \dots, F_{c}^{P}, \dots, F_{m}^{P} {, t]}^{T}

(20)

where $F_{c}^{P}$ is the parameter of fluctuation $c$ ; and the subscript $c \in ℕ^{+}$ indexes each distinct fluctuation among the m-dimensional fluctuations.

When addressing the envelope problem, the proposed fast solution method prioritizes exploring critical regions located on 2-dimensional projections. The proposed fast solution method determines set-points by fixing certain fluctuations and decomposing the m-dimensional $F^{P}$ into $m$ 2-dimensional $F_{c}^{P}$ .

θ = [\begin{array}{l} F^{P} \\ t \end{array}] \Rightarrow θ_{1} = [\begin{matrix} F_{1}^{P} \\ t \end{matrix}], \dots, θ_{c} = [\begin{matrix} F_{c}^{P} \\ t \end{matrix}], \dots, θ_{m} = [\begin{matrix} F_{m}^{P} \\ t \end{matrix}]

(21)

To distinguish among the $m$ 2-dimensional problems, we also use the subscript $c$ to index them. The 2-dimensional parametric programming model is presented as:

\underset{G (θ_{c})}{m i n} z = H_{1} G (θ_{c}) + H_{2} e_{G}^{T}

(22)

s.t.

\{\begin{array}{l} A G (θ_{c}) \leq B_{c} + C_{c} F_{c}^{P} + t D \\ G (θ_{c}) \in 𝒢 \subseteq ℝ^{n}, F_{c}^{P} \in ℱ_{c} \subseteq ℝ^{1}, t \in 𝒯 \subseteq ℝ^{1} \\ θ_{c} : = [F_{c}^{P} {, t]}^{T} θ_{c} \in Θ_{c} \subseteq ℝ^{2} \end{array}

(23)

where $C_{c}$ is a $p \times 1$ constant matrix corresponding to the $c^{t h}$ column of $C$ ; $B_{c}$ is a $p \times 1$ constant matrix that modifies the original $B$ by accounting for other fixed fluctuations (other parameters in $θ$ except $F_{c}^{P}$ ); $ℱ_{c}$ is the 1-dimentional polytope of $F_{c}^{P}$ ; and $Θ_{c}$ is the 2-dimensional polytope of $θ_{c}$ .

By employing the 2-dimensional decompositions in (21), we significantly reduce the computational complexity of the original model in (9). This acceleration is achieved by excluding critical regions that are not present in the 2-dimensional projections. For the m-dimensional $F^{P}$ , the model in (22) and (23) needs to calculate $m$ 2-dimensional projections.

In Fig. 3, the results is visualized using the 2-dimensional fluctuation parameter $F^{P} = [F_{1}^{P}, F_{2}^{P}]^{T}$ and t. In Fig. 3(a), the 2-dimensional critical regions of $F_{2}^{P}$ with t are presented. These critical regions, represented by two different colors, merge to form a 2-dimensional polytope $Θ_{2} (Θ_{c})$ , outlined by the red line. For comparison, the blue dotted line indicates the exact 3-dimensional polytope $Θ$ . The red points mark the vertices of $Θ_{2}$ , while the blue points correspond to the vertices of $Θ$ . In Fig. 3(b), the results of different parameters $F_{1}^{P}$ and $F_{2}^{P}$ are presented using the model in (22) and (23). Two 2-dimensional polytopes, $Θ_{1}$ and $Θ_{2}$ , are obtained, and their vertices lie within $Θ$ . The red dotted line indicates the intersection of the two 2-dimensional projections.

Fig. 3 2-dimensional projections and accurate 3-dimensional polytope. (a) 2-dimensional projection of $Θ_{2}$ . (b) 2-dimensional projections of $Θ_{1}$ and $Θ_{2}$ .

C. Algebraic Manipulation for Original-dimensional Results

Utilizing sensitivity analysis [

28] and parametric programming theory [31], [33], we can efficiently transform the 2-dimensional projections back to the original dimensions without additional optimizations. Algebraic manipulations are used for dimension-raising to obtain the approximated polytope

Θ^{a}

in the original-dimensional parameter space.

In a general parametric programming problem in (9), the Jacobian matrices of the system for decision variables can be derived. These matrices are transformed into matrices $M_{i}$ and $N_{i}$ for the parameter vectors using parametric programming theory [

31].

M_{i}

is an

(n + p) \times (n + p)

matrix and

N_{i}

is an

(n + p) \times (m + 1)

matrix, which represent the Jacobian of the system corresponding to the critical region

ℋ_{i}

. These matrices can be derived with

θ^{*}

as:

\{\begin{array}{l} M_{i} = [\begin{matrix} Q & A_{1}^{T} & \dots & A_{r}^{T} & \dots & A_{p}^{T} \\ - λ_{1} (θ^{*}) A_{1} & - W_{1} & \dots & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ - λ_{r} (θ^{*}) A_{r} & 0 & \dots & - W_{r} & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ - λ_{p} (θ^{*}) A_{p} & 0 & \dots & 0 & \dots & - W_{p} \end{matrix}] \\ N_{i} = [Y^{T}, λ_{1} (θ^{*}) K_{1}^{T}, \dots, λ_{r} (θ^{*}) K_{r}^{T}, \dots, λ_{p} (θ^{*}) K_{p}^{T}]^{T} \\ W_{r} = A_{r} G (θ^{*}) - B_{r} - K_{r} θ^{*} \end{array}

(24)

where Q is an $n \times n$ symmetric constant matrix; Y is an $n \times (m + 1)$ null matrix; $G (θ^{*})$ is the optimal value obtained from the model in (9) at $θ^{*}$ ; $K = [C ∣ D]$ is a $p \times (m + 1)$ matrix, formed by horizontally concatenating the $p \times m$ matrix C and the $p \times 1$ matrix D; the subscript $r \in ℕ^{+}$ refers to the $r^{t h}$ row vector in the original matrices, indexing the $r^{t h}$ constraint among the total p constraints; and $λ_{r} (θ^{*})$ is the Lagrange multiplier of the $r^{t h}$ constraint.

According to (24), once $θ^{*}$ and its corresponding values $G (θ^{*})$ and $λ (θ^{*})$ are obtained, $M_{i}$ and $N_{i}$ for $ℋ_{i}$ can be calculated directly. The algebraic manipulations start by obtaining $θ^{*}$ , $G (θ^{*})$ , and $λ (θ^{*})$ from low-dimensional projections.

To clarify, we add the subscript c to indicate the 2-dimensional results derived from the model in (22) and (23), and the subscript $k \in ℕ^{+}$ indexes subsequent results derived from 2-dimensional critical regions $ℋ_{c, k}$ . We then use algebraic manipulation to transform these critical regions back into their original $ℋ_{k}$ . For each 2-dimensional critical region $ℋ_{c, k}$ , we arbitrarily select a point $θ_{c, k}^{*}$ on $ℋ_{c, k}$ as the set-point. The values of $G (θ_{c, k}^{*})$ and $λ (θ_{c, k}^{*})$ can be directly calculated from the affine functions of parametric solutions associated with $ℋ_{c, k}$ .

Notably, $ℋ_{c, k}$ is a projection of $ℋ_{k}$ . We set values of $θ^{*}$ , $G (θ^{*})$ , and $λ (θ^{*})$ as $θ_{c, k}^{*}$ , $G (θ_{c, k}^{*})$ , and $λ (θ_{c, k}^{*})$ , respectively. Using $θ^{*}$ , $G (θ^{*})$ , and $λ (θ^{*})$ according to (24), we can calculate $M_{k}$ and $N_{k}$ for $ℋ_{k}$ .

After obtaining $M_{k}$ and $N_{k}$ , we calculate the affine functions $G_{k}^{P} (θ)$ and $λ_{k}^{P} (θ)$ associated with $ℋ_{k}$ . As mentioned in Section III-B, in the 2-dimensional results, only a single fluctuation is treated as a varying parameter, while the other fluctuations in the $(m + 1)$ -dimensional vector θ are held constant. When calculating the original-dimensional results, the previously constant fluctuations are treated as varying parameters. $G_{k}^{P} (θ)$ and $λ_{k}^{P} (θ)$ are derived in the neighborhood of $θ^{*}$ as:

[\begin{matrix} G_{k}^{P} (θ) \\ λ_{k}^{P} (θ) \end{matrix}] = - M_{k}^{- 1} N_{k} (θ - θ^{*}) + [\begin{matrix} G (θ^{*}) \\ λ (θ^{*}) \end{matrix}]

(25)

$G_{k}^{P} (θ)$ and $λ_{k}^{P} (θ)$ can further define $ℋ_{k}$ as shown in (12)-(14). Finally, $G_{k}^{P} (θ)$ , $λ_{k}^{P} (θ)$ , and $ℋ_{k}$ are all obtained through algebraic manipulation.

Repeating algebraic manipulations for all 2-dimensional critical regions results in a total of $𝒦$ dimension-raising critical regions. Each $ℋ_{k}$ has a corresponding set of vertices, denoted by $𝒱_{k}$ . These sets are collectively expressed as ${𝒱_{k}}^{𝒦} = {𝒱_{1}, \dots, 𝒱_{k}, \dots, 𝒱_{𝒦}}$ , where $𝒦$ is less than or equal to the total number $ℐ$ of original-dimensional critical regions. Finally, $Θ^{a} : = c o n v ({𝒱_{k}}^{𝒦})$ is defined via V-rep.

Figure 4 builds upon the 2-dimensional projections shown in Fig. 3 to illustrate the transformation into 3-dimensional results.

Fig. 4 Algebraic manipulations for approximated polytope. (a) Dimension-raising operations for 2-dimensional results. (b) Critical regions obtained from dimension-raising and vertices of 3-dimensional approximated polytope.

Notably, some 2-dimensional projections represent different views of the same 3-dimensional critical regions. In Fig. 4(a), all 2-dimensional critical regions are transformed into 3-dimensional critical regions, and consistent colors are used to highlight the corresponding 3-dimensional critical regions. Figure 4(b) removes the 2-dimensional projections to clearly display the 3-dimensional results. The vertices (red points) form $Θ^{a}$ . $Θ^{a}$ , as a subset of the exact polytope $Θ$ , omits certain vertices (blue points) of $Θ$ , indicating accuracy loss when relying solely on low-dimensional projections.

D. Polytope Expansion Strategy for More Accurate Results

We propose a polytope expansion strategy to identify critical regions beyond $Θ^{a}$ . By subdividing the unexplored regions and merging newly identified critical regions, we can achieve a more accurate polytope.

As defined in (14), each $ℋ_{k}$ in the original dimension is determined by a set of halfspaces corresponding to parametric inactive constraints $h_{I, k}$ . At a specific boundary facet of $ℋ_{k}$ , a parametric constraint transitions from inactive to active. This transition is also reflected in the decision variable space, where the corresponding constraint moves from inactive to active. These transitions form the basis for further analysis.

In Section III-C, we have identified several critical regions ${𝒱_{k}}^{𝒦}$ that shape $Θ^{a}$ . Based on ${𝒱_{k}}^{𝒦}$ , we can find boundary facets of $Θ^{a}$ . Each boundary facet of $Θ^{a}$ corresponds to a halfspace. Moving along the normal vector from the interior of $Θ^{a}$ , a parametric constraint transitions from inactive to active upon crossing the boundary facet. Then, at an arbitrary proximal-point just outside this boundary facet, the active and inactive constraints in the decision variable space can be identified based on the corresponding boundary facet, regardless of whether a constraint transition occurs in the decision variable space.

More specifically, consider a parametric constraint j, where $j \in ℕ^{+}$ indexes the boundary facet of $Θ^{a}$ and corresponds to a halfspace of $ℋ_{k}$ . The set of parametric inactive constraints associated with $ℋ_{k}$ , denoted by $h_{I, k}$ , includes the parametric constraint j. Thus, the boundary facet defines a halfspace $𝒫_{k, j}$ . When crossing the boundary facet defined by parametric constraint j, $𝒫_{k, j}$ is flipped, transforming it into the halfspace $𝒬_{k, j}$ .

\{\begin{array}{l} 𝒫_{k, j} : = {θ \in ℝ^{m + 1} : h_{I, k, j} (θ) \leq 0} \\ 𝒬_{k, j} : = {θ \in ℝ^{m + 1} : h_{I, k, j} (θ) \geq 0} \end{array}

(26)

Replacing $𝒫_{k, j}$ with $𝒬_{k, j}$ in $Θ^{a}$ delineates an unexplored region defined by the parametric constraint j. Our goal is to identify the omitted critical regions within this unexplored region to refine the polytope. When the boundary facet defined by parametric constraint j is reached, a previously inactive constraint may become active in the decision variable space, which can be identified by introducing $G_{k}^{P} (θ)$ in the parametric constraint j. At the same time, the remaining constraints can also be confirmed to remain inactive. If all previously inactive constraints remain inactive, the active and inactive constraints in the decision variable space are also clearly identified. Therefore, once the unexplored region is defined, the active and inactive constraints are determined with the introduction of $G_{k}^{P} (θ)$ .

Given that the sets of active and inactive constraints are already established, we now focus on the feasible domain at a given proximal-point. If a feasible domain is confirmed, this proximal-point, located just outside the boundary facet, is selected as the new set-point $θ^{*}$ . $θ^{*}$ is then used to identify critical regions beyond the approximated polytope. This process is repeated for other parametric constraints, defining the unexplored regions and facilitating the discovery of new critical regions.

The advantage of polytope expansion strategy lies in pre-defining unexplored regions in the parameter space and establishing the corresponding sets of active and inactive constraints. For a specific set-point, effective optimization techniques, such as the active-set method [

34] and proximal-point iterations [35], can be employed to navigate the feasible domain. However, as new critical regions emerge, the computational burden of calculating unexplored regions grows, leading to the curse of dimensionality. Although the polytope expansion strategy helps define new regions, it does not inherently resolve the computational challenges associated with high-dimensional problems. Therefore, focusing on specific regions of interest helps strike a balance between polytope accuracy and computational efficiency.

Figure 5 illustrates the polytope expansion strategy for more accurate results, based on $Θ^{a}$ , as shown in Fig. 4. Initially, in Fig. 5(a), $Θ^{a}$ is formed by several halfspaces. Moving along the normal vectors to the boundary facets (indicated by blue arrows), we identify unexplored regions. Among these, only two regions (outlined by red facets) have feasible domains in the decision variable space after flipping all parametric constraints associated with the boundary facets. By optimizing at two specific set-points (proximal points), two new critical regions are discovered within these unexplored regions, as shown in Fig. 5(b). Finally, the polytope expansion strategy allows $Θ^{a}$ to merge these newly identified critical regions, refining the polytope for a more accurate representation of the continuous operating envelope.

Fig. 5 Polytope expansion strategy for more accurate results. (a) Definition of unexplored regions. (b) Identification of new critical regions.

IV. Case Studies

This paper introduces the concept of a continuous operating envelope to support fluctuation management and demonstrates its effectiveness in two power systems: an illustrative 5-bus system and a utility-scale 661-bus system. A 3-dimensional parameter space is visualized to provide an intuitive representation. Furthermore, this paper evaluates the computational performance in handling high-dimensional envelope problems, with a particular focus on the utility-scale 661-bus system results. The findings show that the multi-parametric toolbox (MPT3) [

36] encounters significant difficulties in generating polytopes for higher-dimensional cases, primarily due to the increased dimensionality and its computational burden. In contrast, the proposed fast solution method consistently generates a continuous operating envelope through an approximated polytope, while maintaining predictable computational efficiency.

All optimizations are performed in a MATLAB environment using CPLEX v12.7.1 on a ThinkPad X1 2021 with an Intel^(R) Core^(TM) i5-1135G7 CPU. The benchmark for multi-parametric programming solutions is established using the publicly available MPT3 tool.

A. Illustration of Continuous Operating Envelope

The concept of the continuous operating envelope is detailed in an illustrative 5-bus system, as shown in Fig. 6.

Fig. 6 An illustrative 5-bus system.

This system consists of thermal units (Gen₁ and Gen₂), which are dispatchable resources, and wind units (Gen₃ and Gen₄), which are non-dispatchable resources. The physical limits of the units are listed in Table I. The transmission capacity of Line₄ is 200 MW, while other lines each have a capacity of 400 MW. The reactance of all lines is assumed to be identical.

TABLE I Physical Limits of Units in Illustrative 5-bus System

Resource	Unit	Production cost bid ($/MWh)	Bid capacity (MW)		R_U (MW/h)	R_D (MW/h)
Resource	Unit	Production cost bid ($/MWh)	Maximum	Minimum	R_U (MW/h)	R_D (MW/h)
Gen₁	Thermal	25	700	200	100	100
Gen₂	Thermal	30	500	200	50	50
Gen₃	Wind	0	120	0	-	-
Gen₄	Wind	0	80	0	-	-

Initially, a discrete dispatch decision is made to match the demand over an hourly interval. Specifically, Load₁ and Load₂ are 418.90 MW and 214.82 MW, respectively. Gen₁ and Gen₂ generate 233.72 MW and 200 MW, respectively, while Gen₃ and Gen₄ generate 120 MW and 80 MW, respectively. These values represent the known dispatch results at the start of the interval $(t = t^{S})$ .

To analyze the system ability and handle intra-interval fluctuations, the outputs of Gen₃ and Gen₄ are treated as fluctuating parameters, denoted by $F^{P} = [F_{1}^{P}, F_{2}^{P}]^{T}$ . For simplicity, the loads are assumed to remain constant throughout the intra-interval (intra-hourly) period, without any fluctuations. The goal is to construct a 3-dimensional continuous operating envelope for the output of Gen₃ and Gen₄, spanning from $t^{S} = 0$ min to $t^{E} = 60$ min, based on the known dispatch decision at $t = t^{S}$ .

Assume that Gen₃ and Gen₄ have the installed capacities of 300 MW each, which exceed their maximum bid capacities shown in Table I. This means that Gen₃ and Gen₄ have the potential to increase their generation during the interval. These excess capacities represent potential intra-interval fluctuations due to production uncertainty. In the following analysis, $F_{1}^{P}$ and $F_{2}^{P}$ represent the outputs of Gen₃ and Gen₄, respectively. The difference between the known dispatch result (120 MW for Gen₃ and 80 MW for Gen₄) at the start of the interval $(t = t^{S})$ and the outputs represents deviations within the interval $𝒯$ , defined by $t^{S} \leq t \leq t^{E}$ . For general consideration, $F_{1}^{P}$ and $F_{2}^{P}$ at $t = t^{S}$ are also treated as potential fluctuations. To simplify the subsequent visualization of Gen₃ and Gen₄, we regard the fluctuations $F^{P}$ as outputs of renewables here, without subtracting the constant outputs at $t = t^{S}$ as done in (7).

In Fig. 7(a), the closed lines outline two critical regions in the 2-dimensional space, with each region defined by a distinct set of inactive constraints. Merging these critical regions forms the 2-dimensional polytope $Θ_{1}$ , which represents the continuous operating envelope for Gen₃ in 2-dimensional parameter space. The H-rep of $Θ_{1}$ is given in (27).

Fig. 7 Continuous operating envelope in an illustrative 5-bus system. (a) 2-dimensional projection of Gen₃. (b) 2-dimensional projection of Gen₄. (c) Algebraic manipulation. (d) 3-dimensional polytope.

\{\begin{array}{l} θ_{1} : = [F_{1}^{P} {, t]}^{T} θ_{1} \in Θ_{1} \subseteq ℝ^{2} \\ Θ_{1} : = \{[\begin{matrix} F_{1}^{P} \\ t \end{matrix}] \in ℝ^{2} : [\begin{matrix} \begin{matrix} - 0.4 \\ \begin{matrix} 0.6 \\ 1 \\ 0 \\ - 1 \end{matrix} \end{matrix} & \begin{matrix} - 1 \\ \begin{matrix} - 1 \\ 0 \\ 1 \\ 0 \end{matrix} \end{matrix} \end{matrix}] [\begin{matrix} F_{1}^{P} \\ t \end{matrix}] \leq [\begin{matrix} - 48 \\ \begin{matrix} 72 \\ 153.72 \\ 60 \\ 0 \end{matrix} \end{matrix}]\} \end{array}

(27)

The first and second halfspaces in (27) represent the parametric formulation of ramping constraints, modeled with respect to $t$ . The ramping capacities of Gen₁ and Gen₂ define the feasible operational range at various points in time. The third constraint is a parametric formulation of the transmission in (27), which limits injected power to manage fluctuations. The fourth and fifth halfspaces in (27) represent the parameter bounds. Figure 7(b) presents a similar result for Gen₄, where the 2-dimensional continuous operating envelope, denoted by $Θ_{2}$ , is calculated.

Next, algebraic manipulation is used to extend the previously obtained 2-dimensional projections to the original-dimensional regions. Figure 7(c) demonstrates the process of extending the results to higher dimensions, where the 2-dimensional results are used to construct an approximated polytope. For example, $Θ_{1}$ is extended by calculating the system’s Jacobian corresponding to the critical regions, represented by yellow and blue in Fig. 7(a), respectively. A similar algebraic manipulation is applied to $Θ_{2}$ . Finally, the 3-dimensional envelope is represented by the approximated polytope as shown in Fig. 7(d), with its H-rep given in (28).

\{\begin{array}{l} θ : = [F_{1}^{P}, F_{2}^{P} {, t]}^{T} θ \in Θ \subseteq ℝ^{3} \\ Θ : = \{[\begin{matrix} F_{1}^{P} \\ F_{2}^{P} \\ t \end{matrix}] \in ℝ^{3} : [\begin{matrix} 3 & 3 & - 5 \\ \begin{matrix} - 2 \\ 1 \\ - 1 \\ 0 \end{matrix} & \begin{matrix} - 2 \\ 1 \\ 0 \\ - 1 \end{matrix} & \begin{matrix} - 5 \\ 0 \\ 0 \\ 0 \end{matrix} \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} F_{1}^{P} \\ F_{2}^{P} \\ t \end{matrix}] \leq [\begin{matrix} 600 \\ \begin{matrix} - 400 \\ 233.72 \\ 0 \\ 0 \end{matrix} \\ 60 \end{matrix}]\} \end{array}

(28)

To obtain a more accurate polytope, we continue to use the polytope expansion strategy. Note that the fourth to sixth halfspaces in (28) represent the parameter bounds, which cannot be flipped. The remaining three parametric constraints can be flipped and require additional analysis.

For example, we analyze the first halfspace of $Θ^{a}$ in (28), corresponding to the yellow critical region in Fig. 7(c). The affine functions for Gen₁ $G_{G 1}^{P} (\cdot)$ and Gen₂ $G_{G 2}^{P} (\cdot)$ in this critical region are given as:

G^{P} (θ) = [\begin{matrix} G_{G 1}^{P} (θ) \\ G_{G 2}^{P} (θ) \end{matrix}] = [\begin{matrix} \begin{matrix} - 1 \\ 0 \end{matrix} & \begin{matrix} - 1 \\ 0 \end{matrix} & \begin{matrix} 0 \\ 0 \end{matrix} \end{matrix}] [\begin{matrix} F_{1}^{P} \\ F_{2}^{P} \\ t \end{matrix}] + [\begin{matrix} 433.72 \\ 200 \end{matrix}]

(29)

By substituting the affine functions from (29) into the first halfspace of $Θ^{a}$ , we derive the corresponding constraint expression in the decision variable space as:

233.72 - \frac{5}{3} t \leq G_{G 1}^{P} (θ)

(30)

The constraint (30) represents the ramping constraint, specifically focusing on the ramp-down capacity of Gen₁. At this stage, Gen₂ has already reached its minimum output and cannot ramp down any further. As a result, the first halfspace of $Θ^{a}$ in (28) is flipped to form halfspace $𝒬$ in the parameter space as:

𝒬 = {θ \in ℝ^{3} : 3 F_{1}^{P} + 3 F_{2}^{P} - 5 t \geq 600}

(31)

We select a proximal-point (set-point) of $F_{1}^{P} = 120$ MW, $F_{2}^{P} = 101$ MW, and $t = 12$ min as an example. This point satisfies (31) and is located outside the boundary facet of (28). At $t = 12$ min, it shows that $G_{G 1}^{P} (θ)$ can only decrease by 20 MW, which is insufficient to accommodate the 21 MW increase in Gen₄ output. In fact, we can validate that the system is not feasible within the halfspace defined by (31), as the inequality constraint in (30) becomes active, and the ramping-down capability of the system is exhausted. Therefore, no feasible domain exists after flipping the boundary facet indicated by the first halfspace of $Θ^{a}$ in (28).

When performing the same operation on the second and third halfspaces of $Θ^{a}$ in (28), we identify no unexplored regions in the parameter space. This is due to the exhausted ramping and transmission capacities in the decision space. As a result, no further optimization or expansion of the polytope is needed. The polytope $Θ^{a}$ in Fig. 7(d) accurately represents the continuous operating envelope. Here, the approximated polytope $Θ^{a}$ is accurate due to the small number of parametric constraints. A detailed discussion of the approximated polytope is provided in the subsequent case. Moreover, a direct solution using MPT3 yields the same polytope, as shown in Fig. 7(d).

In summary, the proposed fast solution method incorporates the fluctuations of two wind units and t to determine a 3-dimensional polytope, which represents the continuous operating envelope. This polytope is shaped by constraints in (9), with key factors being the ramping capacities of dispatchable resources (Gen₁ and Gen₂) and the transmission capacity.

B. Efficiency of Continuous Operating Envelope

To demonstrate the efficiency of the continuous operating envelope in managing fluctuations, we begin by analyzing the intra-interval fluctuations of Gen₃ and Gen₄, each represented by its respective 2-dimensional envelope. Assuming that fluctuations are recorded every 15 s, this results in 240 set-points (as shown in Fig. 8), forming a fluctuation curve over an hour. These data are based on actual wind unit samples. As shown in Fig. 8(a) and Fig. 8(b), all individual fluctuations fall within their respective 2-dimensional envelopes, indicating that the fluctuations are manageable.

Fig. 8 Fluctuation management using continuous operating envelope. (a) Validation of Gen₃ fluctuations. (b) Validation of Gen₄ fluctuations. (c) Validation of simultaneous fluctuations of Gen₃ and Gen₄. (d) Identification of infeasible points among nine additional curves of simultaneous fluctuations.

However, when Gen₃ and Gen₄ fluctuate simultaneously, the interaction between their fluctuations leads to intra-interval infeasibility. Consequently, some set-points become infeasible. As shown in Fig. 8(c), the envelope defined in (28) clearly distinguishes between feasible set-points (marked in green) and infeasible ones (marked in red). For comparison, the 3-dimensional infeasible fluctuations are also projected as red points in Fig. 8(a) and Fig. 8(b).

The continuous operating envelope is particularly effective in directly pinpointing infeasible fluctuations. In the scenarios where Gen₃ and Gen₄ frequently adjust their outputs in real time, feasibility examinations become increasingly complex. Each update involves checking the feasibility of 240 set-points, which becomes cumbersome when performed repeatedly. For example, examining nine fluctuation curves would require performing 240 optimizations for each curve. In contrast, the continuous operating envelope simplifies this process by directly identifying infeasible fluctuations (marked in red), as demonstrated in Fig. 8(d).

Thus, the continuous operating envelope not only streamlines the assessment of operational risks due to intra-interval fluctuations but also significantly reduces the workload involved in performing repeated feasibility examinations for each update.

C. Computational Performance of Proposed Fast Solution Method

We assess the computational performance of the proposed fast solution method on a utility-scale 661-bus system, which includes 37 generation units and 1047 transmission lines. In this paper, we focus initially on two wind units as primary sources of intra-interval fluctuations and calculate their corresponding continuous operating envelope.

In Fig. 9, the results generated by the proposed fast solution method are illustrated. As depicted in Fig. 9(a) and Fig. 9(b), two 2-dimensional projections are calculated, with the black regions indicating the accumulation of critical regions. These 2-dimensional projections are then algebraically transformed into their original 3-dimensional regions. In different 2-dimensional projections, critical regions that correspond to the same 3-dimensional regions are shown in the same color. For clarity, Fig. 9(c) highlights two representative critical regions in purple and blue, demonstrating how the dimension-raising process works. The final approximated polytope is outlined in blue in Fig. 9(c) after the algebraic manipulation process.

Fig. 9 Continuous operating envelope in utility-scale 661-bus system. (a) 2-dimensional projection of Wind₁. (b) 2-dimensional projection of Wind₂. (c) 3-dimensional polytope.

However, as the number of parametric constraints increases in more complex systems, such as the utility-scale 661-bus system, the approximated polytope obtained from low-dimensional projections may result in accuracy loss. To mitigate this, the polytope expansion strategy can be applied to enhance accuracy. By flipping the boundary facets of the approximated polytope (outlined in blue), unexplored regions in the parameter space are revealed. In Fig. 9(c), two such unexplored regions (in red) are outlined, where feasible domains exist. To refine the polytope and achieve higher accuracy, these unexplored regions need to be further subdivided and explored.

In particular, the unexplored region in the upper left corner of Fig. 9(c) is selected for detailed analysis. As shown in Fig. 10(a), this region is subdivided by calculating a sufficient number of set-points to explore its full extent. Through the polytope expansion strategy, these subdivisions reveal new critical regions beyond the original approximated polytope. Once all unexplored regions are fully subdivided, the accurate polytope is constructed by merging the new critical regions with the initial approximated polytope. This refined polytope, as outlined in blue in Fig. 10(b), more accurately represents the continuous operating envelope.

Fig. 10 Continuous operating envelope obtained by proposed fast solution method. (a) Illustration of polytope expansion strategy. (b) Accurate polytope.

Table II highlights the computational efficiency for 3-dimensional results. In the 2-dimensional projections of different fluctuations, some critical regions are the projections of the same 3-dimensional critical regions. Here, 48 critical regions are shared.

TABLE Ⅱ Computational Efficiency for 3-dimensional Results

Method	Number of explored critical regions	Volume of polytope (%)	Time (s)
2-dimensional projections (Wind₁)	49	-	29.12
2-dimensional projections (Wind₂)	61	-	36.26
Algebraic manipulation	62	99.42	0
Polytope expansion strategy	553	100.00	638.65
Benchmark (MPT3)	659	100.00	734.83

Besides, the time required for algebraic manipulation to generate the matrices is less than 0.01 s. The approximated polytope (algebraic manipulation part) is constructed with 62 explored critical regions, which results in a slight reduction in volume and indicates a loss in accuracy. To further refine the results, the polytope expansion strategy is applied, leading to the identification of the exact polytope with 553 explored critical regions. Compared with the benchmark, the proposed fast solution method efficiently constructs the exact polytope, excluding 106 non-outermost critical regions.

As shown in Table II, the polytope expansion strategy is the most time-consuming, which takes 638.65 s. Despite the extended time, the volume of the polytope increases by only 0.58%, indicating that many critical regions in the unexplored region contribute minimally to the overall volume but significantly increase the computational burden. In contrast, an approximated polytope can be computed in just 65.38 s before applying the polytope expansion strategy. It shows that the polytope expansion refines accuracy without significantly altering the overall polytope.

With sufficient time, the proposed fast solution method can generate a comprehensive and exact polytope. However, the computational burden of the polytope expansion strategy remains on the same order of magnitude as MPT3 (benchmark). In time-constrained scenarios, the proposed fast solution method provides a practical trade-off by generating an approximated polytope, either without or with partial polytope expansion. While this method may incur some loss in accuracy, it significantly reduces the computational burden.

To further validate the accuracy of the envelope, we incorporate four additional wind units at different buses. Each new wind unit introduces a parametric fluctuation, increasing the dimensionality of the continuous operating envelope (six parametric fluctuations in total, plus one parametric timescale). In this analysis, low-dimensional projections are used to handle high-dimensional results. Previously calculated low-dimensional information remains valid when new fluctuations are introduced, with only the newly fluctuating nodes requiring recalculation. This new information is then algebraically transformed into the high-dimensional polytope.

For example, when a new fluctuating node is introduced into a 2-dimensional parameter space, the dimensionality increases to 3. As shown in Table III, the additional time required is $65.38 - 29.12 = 36.26$ s, which corresponds to calculating the 2-dimensional projection of the new fluctuating node. Table III also shows that the total parameter dimensions reach 7 by calculating six independent 2-dimensional projections. The total computation time is the sum of these independent 2-dimensional results. If previous results are available, updating the continuous operating envelope (approximated polytope) only requires calculating the additional 2-dimensional projections of newly fluctuating nodes. Moreover, parallel algorithms can be used to speed up the process when multiple fluctuations occur simultaneously.

TABLE Ⅲ Performance of Approximated Polytope and MPT3 in Various Dimensions

Dimension	Approximated polytope			MPT3
Dimension	Number of critical regions	Time (s)	Volume (%)	Number of critical regions	Time (s)
2	49	29.12	100.00	49	29.12
3	62	65.38	99.42	659	734.83
4	78	82.34	96.99	1456	1386.47
5	81	93.47	98.64	1682	1683.62
6	108	119.06	-	-	-
7	121	144.11	-	-	-

In high-dimensional problems, where volume is difficult to intuitively define, envelope accuracy is estimated by comparing set-points to a benchmark. For instance, we use the 4-dimensional MPT3 results in Table III as the benchmark, generating 5 million set-points within the accurate polytope to verify them against the approximated polytope. This sampling method is also employed to estimate the volume discrepancy of the 5-dimensional results. Note that in Table III, the 4-dimensional sampling result is 4849340/5000000; and the 5-dimensional sampling result is 4931986/5000000. The volume of the approximated polytope closely matches that of the accurate MPT3 polytope, while significantly reducing the computational burden.

Table III highlights the performance of the approximated polytope across different dimensions. The results show that computational cost increases linearly with dimensionality, while accuracy loss remains limited. This demonstrates that the proposed fast solution method can stably generate an approximated polytope, offering a significant advantage over MPT3 in higher-dimensional cases. For instance, in the 6-dimensional case, after evaluating 3500 critical regions, the number of regions has yet converged, leading us to conclude that the problem is unsolvable using MPT3.

The discussion above focuses on the approximated polytope without the polytope expansion strategy. In time-constrained scenarios, the approximated polytope can still serve as a continuous operating envelope, with fluctuations within the envelope considered feasible. Fluctuations outside the envelope, while potentially feasible, are treated as operational risks due to the exclusion of high-dimensional critical regions, which results in some accuracy loss. Although the polytope expansion strategy can improve accuracy, it is time-consuming for high-dimensional problems, much like traditional parametric programming methods. When time permits, the polytope expansion strategy can be employed to obtain a more precise polytope, but its computational burden increases exponentially with dimensionality.

In summary, in time-constrained situations, the approximated polytope offers a rapid solution that effectively represents the continuous operating envelope for intra-interval fluctuations. Conversely, when more time is available, the polytope expansion strategy can be applied to obtain a more accurate polytope. This method is particularly advantageous in high-dimensional problems, ensuring that at least an approximated polytope is attained to conservatively manage fluctuations. In practice, when frequent updates to the continuous operating envelope are necessary, a conservative approximated polytope can serve as the envelope, with only low-dimensional analysis needed for newly added fluctuating nodes.

V. Conclusion

In power systems with discrete interval horizons, maintaining intra-interval feasibility is a challenging task due to fluctuations. This paper introduces the concept of a continuous operating envelope, which provides physical insights into the allowable range of intra-interval fluctuations. The continuous operating envelope is formulated as a parametric programming model, and its computational complexity is theoretically analyzed. To overcome the inherent computational challenges, we develop a fast solution method to efficiently construct the envelope and ensure the scalability for high-dimensional problems. The proposed fast solution method is validated through both an illustrative 5-bus system to illustrate the concept and a utility-scale 661-bus system to verify its scalability and efficiency.

The continuous operating envelope offers a proactive tool for fluctuation management, enabling more reliable and efficient power system operations. By defining a clear range of manageable fluctuations, the proposed fast solution method has the potential to support real-time dispatch decisions and reduce operational risks associated with intra-interval fluctuations.

References

D. S. Kirschen, “Risk exposure, externalization and allocation in unbundled power systems,” IEEE Transactions on Power Systems, vol. 36, no. 3, pp. 1879-1889, May 2021. [Baidu Scholar]

E. Litvinov, F. Zhao, and T. Zheng, “Electricity markets in the United States: power industry restructuring processes for the present and future,” IEEE Power Energy Magazine, vol. 17, no. 1, pp. 32-42, Jan. 2019. [Baidu Scholar]

D. Xiao, H. Chen, W. Cai et al., “Integrated risk measurement and control for stochastic energy trading of a wind storage system in electricity markets,” Protection and Control of Modern Power Systems, vol. 8, no. 4, pp. 1-11, Oct. 2023. [Baidu Scholar]

C. Wu, G. Hug, and S. Kar, “Risk-limiting economic dispatch for electricity markets with flexible ramping products,” IEEE Transactions on Power Systems, vol. 31, no. 3, pp. 1990-2003, May 2016. [Baidu Scholar]

ERCOT. (2012, Jan.). ERCOT fast-responding regulation service, pilot project. [Online]. Available: http://www.ercot.com/mktrules/pilots/frrs [Baidu Scholar]

California ISO. (2017, Jun.). SIBR interface specification for bidding services. [Online]. Available: https://www.caiso.com/Documents/SIBR InterfaceSpecificationWebServicesFall2017Release.pdf [Baidu Scholar]

J. Zhao, T. Zheng, and E. Litvinov, “Variable resource dispatch through do-not-exceed limit,” IEEE Transactions on Power Systems, vol. 30, no. 2, pp. 820-828, Mar. 2015. [Baidu Scholar]

M. Kazemi, P. Siano, D. Sarno et al., “Evaluating the impact of sub-hourly unit commitment method on spinning reserve in presence of intermittent generators,” Energy, vol. 113, pp. 338-354, Oct. 2016. [Baidu Scholar]

M. Parvania and A. Scaglione, “Unit commitment with continuous-time generation and ramping trajectory models,” IEEE Transactions on Power Systems, vol. 31, no. 4, pp. 3169-3178, Jul. 2016. [Baidu Scholar]

M. Zhang, Z. Yang, W. Lin et al., “Enhancing economics of power systems through fast unit commitment with high time resolution,” Applied Energy, vol. 281, p. 116051, Jan. 2021. [Baidu Scholar]

J. Zou, S. Ahmed, and X. Sun, “Multistage stochastic unit commitment using stochastic dual dynamic integer programming,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 1814-1823, May 2019. [Baidu Scholar]

Y. Cho, T. Ishizaki, N. Ramdani et al., “Box-based temporal decomposition of multi-period economic dispatch for two-stage robust unit commitment,” IEEE Transactions on Power Systems, vol. 34, no. 4, pp. 3109-3118, Jul. 2019. [Baidu Scholar]

H. Pandzic, Y. Dvorkin, T. Qiu et al., “Toward cost-efficient and reliable unit commitment under uncertainty,” IEEE Transactions on Power Systems, vol. 31, no. 2, pp. 970-982, Mar. 2016. [Baidu Scholar]

Y. Yu, Security Regions of Power System. Beijing: China Electric Power Press, 2014. [Baidu Scholar]

C. Qin and Y. Yu. “Small signal stability region of power systems with DFIG in injection space,” Journal of Modern Power Systems and Clean Energy, vol. 1, no. 2, pp. 127-133, Sept. 2013. [Baidu Scholar]

F. Wu and Y. Tsai, “Probabilistic dynamic security assessment of power systems: part I-basic model,” IEEE Transactions on Circuits Systems, vol. 30, pp. 148-159, Mar. 1983. [Baidu Scholar]

W. Dai, Z. Yang, J. Yu et al., “Security region of renewable energy integration: characterization and flexibility,” Energy, vol. 187, p. 115975, Nov. 2019. [Baidu Scholar]

T. Yang and Y. Yu, “Static voltage security region-based coordinated voltage control in smart distribution grids,” IEEE Transactions on Smart Grid, vol. 9, no. 6, pp. 5494-5502, Nov. 2018. [Baidu Scholar]

X. Jiang, Y. Zhou, W. Ming et al., “Feasible operation region of an electricity distribution network,” Applied Energy, vol. 331, p. 120419, Feb. 2023. [Baidu Scholar]

Y. Yu and C. Qin, “Security region based security-constrained unit commitment,” Science China Technological Sciences, vol. 56, pp. 2732-2744, Oct. 2013. [Baidu Scholar]

Y. Yu, Y. Liu, C. Qin et al., “Theory and method of power system integrated security region irrelevant to operation states: an introduction,” Engineering, vol. 6, no. 7, pp. 754-777, Jul. 2020. [Baidu Scholar]

S. Vaishya, “An AC optimal power flow framework for active-reactive power scheduling considering generator capability curve,” Energy Conversion and Economics, vol. 4, no. 6, pp. 425-438, Dec. 2023. [Baidu Scholar]

F. Li, T. Niu, L. Xue et al., “Autonomous-synergic voltage security regions in bulk power systems,” Journal of Modern Power Systems and Clean Energy, vol. 11, no. 2, pp. 686-692, Mar. 2023. [Baidu Scholar]

S. Maihemuti, W. Wang, H. Wang et al., “Voltage security operation region calculation based on improved particle swarm optimization and recursive least square hybrid algorithm,” Journal of Modern Power Systems and Clean Energy, vol. 9, no. 1, pp. 138-147, Jan. 2023. [Baidu Scholar]

M. Al-Abdullah, M. Abdi-Khorsand, and W. Hedman, “The role of out-of-market corrections in day-ahead scheduling,” IEEE Transactions on Power Systems, vol. 30, no. 4, pp. 1937-1946, Jul. 2015. [Baidu Scholar]

B. Zad, J. Toubeau, B. Vatandoust et al., “Enhanced integration of flow-based market coupling in short-term adequacy assessment,” Electric Power Systems Research, vol. 201, p. 107507, Dec. 2021. [Baidu Scholar]

CAISO. (2011, Jul.). Dynamic ramp rate in ancillary service procurement. [Online]. Available: https://www.caiso.com/Documents/TechnicalBulletin-DynamicRampRate_AncillaryServiceProcurement.pdf [Baidu Scholar]

A. Fiacco and Y. Ishizuka, “Sensitivity and stability analysis for nonlinear programming,” Annals of Operations Research, vol. 27, pp. 215-236, Dec. 1990. [Baidu Scholar]

J. Ryu, V. Dua, and E. Pistikopoulos, “Proactive scheduling under uncertainty: a parametric optimization approach,” Industrial & Engineering Chemistry Research, vol. 46, no. 24, pp. 8044-8049, Nov. 2007. [Baidu Scholar]

V. Dua, N. Bozinis, and E. Pistikopoulos, “A new multiparametric mixed-integer quadratic programming algorithm,” Computer Aided Chemical Engineering, vol. 9, pp. 979-984, May 2001. [Baidu Scholar]

V. Dua, N. A. Bozinis, and E. N. Pistikopoulos, “A multiparametric programming approach for mixed-integer quadratic engineering problems,” Computers & Chemical Engineering, vol. 26, no. 4-5, pp. 715-733, May 2002. [Baidu Scholar]

P. Orlik and H. Terao, Arrangements of Hyperplanes. Berlin-Heidelberg: Springer-Verlag, 1992. [Baidu Scholar]

V. Sakizlis, K. Kouramas, and E. Pistikopoulos, Linear Model Predictive Control via Multiparametric Programming. New Jersey: Wiley-VCH, 2007. [Baidu Scholar]

P. Hungerländer and F. Rendl, “A feasible active set method for strictly convex quadratic problems with simple bounds,” SIAM Journal on Optimization, vol. 25, no. 3, pp. 1633-1659, Jan. 2015. [Baidu Scholar]

D. Arnstrom, A. Bemporad, and D. Axehill, “Complexity certification of proximal-point methods for numerically stable quadratic programming,” IEEE Control Systems Letters, vol. 5, no. 4, pp. 1381-1386, Oct. 2021. [Baidu Scholar]

M. Herceg, M. Kvasnica, C. Jones et al., “Multi-parametric toolbox 3.0,” in Proceedings of 2013 European Control Conference (ECC), Zurich, Switzerland, Jun. 2013, pp. 502-510. [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