Dispatch Instruction Disaggregation for Virtual Power Plants Using Multi-Parametric Programming

Abstract

Virtual power plants (VPPs) coordinate distributed energy resources (DERs) to collectively meet grid dispatch instructions. When a dispatch command is issued to a VPP, it must be disaggregated optimally among the individual DERs to minimize overall operational costs. However, existing methods for VPP dispatch instruction disaggregation often require solving complex optimization problems for each instruction, posing challenges for real-time applications. To address this issue, we propose a multi-parametric programming-based method that yields an explicit mapping from any given dispatch instruction to an optimal DER-level deployment strategy. In our approach, a parametric optimization model is formulated to minimize the dispatch cost subject to DER operational constraints. By applying Karush–Kuhn–Tucker (KKT) conditions and recursively partitioning the DERs’ adjustable capacity space into critical regions, we derive analytical expressions that directly map dispatch instructions to their corresponding resource allocation strategies and optimal scheduling costs. This explicit solution eliminates the need to repeatedly solve the optimization problem for each new instruction, enabling fast real-time dispatch decisions. Case study results verify that the proposed method effectively achieves the cost-efficient and computationally efficient disaggregation of dispatch signals in a VPP, thereby improving its operational performance.

1. Introduction

In the development of modernized power systems [1,2], the dispersed integration of vast numbers of distributed resources, including distributed generations (DGs), energy storage systems (ESSs), and controllable loads, into distribution systems poses significant operational challenges [3]. As a novel energy management paradigm, virtual power plants (VPPs) aggregate geographically dispersed resources across different grid areas through information and communication technologies (ICT), enabling them to respond to grid dispatching demands as a unified entity. This effectively mitigates the aforementioned pressures [4]. Consequently, coordinating the relationships among resource owners, the VPP, and the distribution network has become a critical aspect of VPP development [5]. However, disaggregating upper-level power commands among heterogeneous internal resources in real time remains a major challenge [6].

In the existing research and practical engineering applications, VPPs primarily employ a direct resource control approach [7] for managing distributed energy resources (DERs). This direct control mode requires resource owners to relinquish their dispatch control rights and benefit allocation rights directly to the VPP [8,9]. The VPP then interfaces directly with the upper-level grid, resulting in a simplified dispatch structure. Under this control paradigm, the VPP operator occupies the core position within the entire dispatch hierarchy. It possesses comprehensive network information, along with the operational status, regulation capabilities, and cost parameters of all subordinate DERs [10]. Reference [11] formulates VPP control as a linear programming problem, designing a real-time energy management system (EMS) for the VPP. This system, coupled with a bidirectional communication infrastructure, enables the optimized control of wind turbines, ESS, and flexible controllable loads. Reference [12] developed a VPP optimal dispatch controller. This controller employs a binary backtracking search algorithm to determine the on/off sequencing of internal units (including PV, wind, ESS, and gas turbines), achieving operational cost minimization. Reference [13] established a dispatch model based on multi-stage distributionally robust optimization (DRO) to decompose the power requirements from the upper-level system into setpoints for individual resources. While this model enhances the dispatch’s adaptability to uncertainties, its architecture mandates that the VPP holds control rights over all resources, and the computational complexity at the central center grows exponentially with the number of resources. Reference [14] aggregates ESS, flexible loads, electric vehicles (EVs), and distributed renewable energy sources into a single VPP entity, which then participates in the grid’s unified dispatch as a whole. However, unified dispatch instructions may induce local network overloading.

Due to high training costs, safety constraint handling errors, and limited generalization capabilities in decentralized AI scheduling frameworks [15], decentralized disaggregating dispatch instructions remain challenging. For reliability, VPPs typically rely on their control centers to calculate disaggregation schedules in real time and execute them directly to DERs. Reference [16] presents a multi-time scale economic scheduling framework for VPP that utilizes the disaggregation of dispatch instructions to enhance flexibility and reduce costs through the intelligent aggregation and disaggregation of deferrable loads. Reference [17] proposes a feature point set-based aggregation method that achieves the efficient disaggregation of dispatch instructions by optimizing the convex combination space of feature point sets, thereby decomposing the aggregate dispatch command received by a VPP into executable sub-instructions for DERs. Reference [18] proposes a coordination operation framework for VPP and distribution network operators. The core innovation lies in achieving the efficient and economical decomposition of scheduling instructions through an optimal disaggregation method based on aggregated cost functions.

For real-time control disaggregation strategies in VPPs, Chen et al. proposed a real-time operational control strategy for VPPs, wherein the power adjustment requested by the grid is disaggregated across heterogeneous internal resources [6]. Liu et al. introduced a multi-time-scale economic dispatch strategy that incorporates the aggregation and disaggregation of deferrable loads into VPP scheduling [18]. This strategy decomposes the grid power demand into different types of controllable resources and coordinates optimization at both day-ahead and real-time scales, thereby enhancing the VPP’s ability to track upper-level commands. Building on this foundation, Huang et al. conducted a comprehensive review of VPP response capability assessment and dispatch optimization techniques [19]. The review scrutinizes the core technologies in VPP dispatch optimization, especially the assessment of aggregated resource response capabilities and optimal allocation methods. It emphasizes that accurately evaluating the controllability of aggregated DERs and optimizing internal resource allocation are key to meeting grid instructions while safeguarding the interests of individual DER units, thus underscoring the significance of the present study. The abovementioned VPP control approaches place a heavy reliance on a central computing node, which requires all operational data to be aggregated for computation and imposes a significant burden on this central node. Particularly in the context of massive DER integration, VPPs necessitate the processing of high-dimensional complex variables and constraints, presenting challenges for the rapid and efficient disaggregation of dispatch instructions to individual resources.

To address the aforementioned challenges, this paper proposes a multi-parametric programming-based method for disaggregating VPP dispatch instructions. Multi-parametric linear programming (MP-LP), as a tool for solving parameterized optimization problems, has made great progress over the past two decades. A hybrid active set strategy improved exploration efficiency [20]. Degeneracy handling enhanced algorithmic robustness [21]. Geometric partitioning via Delaunay triangulation expanded the size of solvable problems [22]. In addition, multi-parametric algorithms that involve mixed-integer or multi-level decision structures, as well as machine learning-assisted approximations, are emerging, offering a solid algorithmic foundation for this study in the use of multi-parametric programming. Applications of MP-LP in power systems are gradually unfolding. In optimal control, explicit MPC and related methods provide analytical control strategies for fast regulation in building energy systems [23]. In uncertainty analysis and planning decisions, MP-LP helps derive functional relationships between input parameters and decisions/costs, greatly improving analytical efficiency [24]. Together, these applications demonstrate that introducing MP-LP into power system optimization can strike a balance between optimality and real-time performance, providing new perspectives for addressing dispatch and control challenges in complex energy systems. In this paper, MP-LP is employed to achieve a real-time mapping from VPP commands to power outputs.

This paper proposes a multi-parametric programming-based method for disaggregating VPP dispatch instructions. The main contributions are as follows:

To resolve the real-time computational bottleneck in scheduling massive distributed resources, a comprehensive strategy and a multi-parametric programming model for VPP dispatch instruction disaggregation are proposed, with the objective of minimizing the total dispatch cost while respecting DER constraints.

By solving the KKT conditions for each critical region, we derive explicit mapping relationships between the VPP’s dispatch instruction and the corresponding resource deployment strategy, as well as the optimal scheduling cost. This provides a closed-form relationship that directly links any dispatch instruction to its optimal disaggregation outcome, ensuring the rapid determination of DER setpoints without the need for iterative optimization.

The rest of this paper is organized as follows: Section 2 presents the framework for VPP dispatch instruction disaggregation. Section 3 presents the multi-parametric programming-based model for VPP dispatch instruction disaggregation. Section 4 develops the computational framework for multi-parametric program transformations and solutions. Section 5 tests the proposed method in case studies. Section 6 provides a summary and concluding remarks.

2. Framework for VPP Dispatch Instruction Disaggregation

The core idea of our method is to shift the solving pressure to the earlier stage. During the day-ahead stage, when time is relatively abundant. The real-time dispatch instructions from the upper grid are treated as parameters. An explicit mapping between the instructions and the optimal dispatch strategy is prepared in advance, allowing the real-time scheduling stage to simply look up the table and execute directly. This process involves techniques like multi-parametric programming based on KKT conditions. Specifically, during the day-ahead stage, (a) active and inactive constraints are differentiated; (b) the critical regions of the parameter domain are recursively partitioned; and (c) a parameterized strategy repository is formed, creating an explicit mapping of the instruction → critical region index → analytical linear strategy/cost.

In the real-time scheduling stage, when the VPP receives the specific “upper grid dispatch instruction,” it looks up the table to identify the corresponding critical region, and based on the explicit mapping of that region, the optimal strategy is determined.

The optimal solution and cost within each critical region are linear functions of the parameters. Therefore, compared to traditional methods that require solving the full optimization model online, our method in the real-time phase only requires O(log N) time (using interval or segment binary search/hash lookup) or O(1) time (direct interval location) to obtain the solution, which is a significant advantage.

This study focuses on VPPs participating in ancillary service markets, shown in Figure 1, which have access to grid information. During the day-ahead stage, the VPP submits a 96-period (15-min interval) power curve (without scheduling), along with upper/lower power output limits for the entire operating day. Leveraging the operational characteristics of its aggregated resources—including flexible loads, DGs, and ESSs—the VPP determines the dispatch schedule associated with its available regulation capacity.

During the real-time stage, the VPP must provide power generation or consumption reduction service, achievable within 10 min, in response to grid dispatch instructions. Failure to execute dispatch commands during real-time operations will result in penalties.

3. Multi-Parametric Programming-Based Model for VPP Dispatch Instruction Disaggregation

3.1. Objective Function

The objective of the VPP instruction disaggregation is to minimize its total operational cost in a single time period, t.

(1)$\min F_{\mathrm{VPP}, t}=C_{\mathrm{GT}, t}+C_{\mathrm{PV}, t}+C_{\mathrm{WT}, t}+C_{\mathrm{ESS}, t}+C_{\mathrm{IL}, t}+C_{\mathrm{penalty}}(P_{\mathrm{Instru} t}-P_{\mathrm{net}, t})$

Here, $F_{\mathrm{VPP}, t}$ denotes the operational cost of the VPP; $C_{\mathrm{GT}, t}$ represents the operational cost of micro gas turbines; $C_{\mathrm{PV}, t}$ is the regulation cost for PV generation; $C_{\mathrm{WT}, t}$ indicates the regulation cost for wind power generation; $C_{\mathrm{ESS}, t}$ is the regulation cost for ESS; $C_{\mathrm{IL}, t}$ defines the compensation cost for demand response loads; and $C_{\mathrm{penalty}}(P_{\mathrm{Instru} t}-P_{\mathrm{net}, t})$ denotes the penalty term of for failing to execute real-time dispatch commands. An additional explanation is as follows:

(2)$C_{\mathrm{penalty}}(P_{\mathrm{Instru}, t}-P_{\mathrm{net}, t})=c_{\mathrm{penalty}}\left|P_{\mathrm{Instru} t}-P_{\mathrm{net}, t}\right|$

where $c_{\mathrm{penalty}}$ denotes the unit price of penalty, which is much bigger than the operational cost, $P_{\mathrm{Instru} t}$ denotes the instruction of the upper-level grid, and $P_{\mathrm{net}, t}$ denotes the actual net power of the distribution system under the disaggregation scheme. The components of the objective function (1) and (2) are defined as follows:

Operational Cost of Micro Gas Turbines

The micro gas turbines managed by the VPP are typically small-capacity units. Considering that operating startup events require additional costs, their operational cost can be modeled as follows:

(3)$\begin{array}{c}{C}_{\mathrm{GT},i,t}=a_i{P}_{\mathrm{GT},i,t}+c_i,P_{\mathrm{GT},i,t-1}>0\\ C_{\mathrm{GT},i,t}=a_i{P}_{\mathrm{GT},i,t}+c_i+C_{\mathrm{GT},i}^{\mathrm{start}},P_{\mathrm{GT},i,t-1}=0\\ C_{\mathrm{GT}, t}={\displaystyle \sum _i{C}_{\mathrm{GT},i,t}}\end{array}$

where $P_{\mathrm{GT},i,t}$ denotes the active power output of the i-th micro gas turbine, $a, c$ correspond to the operational cost coefficients of the gas turbine unit, $P_{\mathrm{GT},i,t-1}$ denotes the deployed scheduling scheme at time interval t-1, and $C_{\mathrm{GT},i}^{\mathrm{start}}$ denotes the additional cost of startup events.

It is worth noting that in the day-ahead calculation, all micro gas turbines are assumed to be off. During intra-day real-time scheduling, if the unit was already started in the instruction decomposition of interval t − 1, the strategy library for the current interval should be recalculated with the first expression of the abovementioned formula before the instruction at time, t, is issued. Because only a constant term is different, the recomputation scale is fully controllable.

2.. Regulation cost of photovoltaic (PV) units

(4)$C_{\mathrm{PV},t}=f_{\mathrm{PV}}\left(P_{\mathrm{PV},t}^{\mathrm{pre}}-P_{\mathrm{PV},t}\right)$

Here, $P_{\mathrm{PV},t}$ denotes the active power output of the PV unit, $f_{\mathrm{PV}}$ corresponds to the operation and maintenance (O&M) cost coefficient of the PV unit, and $P_{\mathrm{PV},t}^{\mathrm{pre}}$ represents the forecast PV power.

3.. Regulation cost of wind turbine units

(5)$C_{\mathrm{WT},t}=f_{\mathrm{WT}}\left(P_{\mathrm{WT},t}^{\mathrm{pre}}-P_{\mathrm{WT},t}\right)$

Here, $P_{\mathrm{WT},t}$ denotes the active power output of the wind turbine unit, $f_{\mathrm{WT}}$ represents the O&M cost coefficient for wind power generation, and $P_{\mathrm{WT},t}^{\mathrm{pre}}$ corresponds to the forecast wind power.

4.. Energy storage regulation cost

(6)$C_{\mathrm{ESS},t}=f_{\mathrm{ESS}}\left|P_{\mathrm{ESS},i,t}\right|$

Here, $C_{\mathrm{ESS}}$ denotes the regulation cost for energy storage, $f_{\mathrm{ESS}}$ denotes the coefficient of regulation cost for energy storage units, and $P_{\mathrm{ESS},i,t}$ denotes the charging and discharging power of the i-th energy storage device during period t.

5.. Demand response load compensation cost

The compensation cost for interruptible loads is formulated as follows:

(7)$C_{\mathrm{IL},t}={\displaystyle \sum _i{c}_{\mathrm{IL}}{P}_{\mathrm{IL},i,t}}$

where $C_{\mathrm{IL}}$ denotes the load reduction compensation expenditure, $c_{\mathrm{IL}}$ denotes the power-based incentive rate, and $P_{\mathrm{IL},i,t}$ denotes the metered curtailable power of participant i during period t.

3.2. Constraint Conditions

Power Balance Constraint

The VPP is required to maintain system-wide power balance at every scheduling period, t. Therefore, the net load of the system under the disaggregation scheme can be expressed mathematically as follows.

(8)$P_{\mathrm{net},t}=P_{\mathrm{L},t}-{\displaystyle \sum _i{P}_{\mathrm{IL},i,t}}-P_{\mathrm{PV},t}-P_{\mathrm{WT},t}-{\displaystyle \sum _i{P}_{\mathrm{GT},i,t}}+{\displaystyle \sum _i{P}_{\mathrm{ESS},i,t}}$

2.. Power output constraints for micro gas turbine units

Given the rapid power response capability of micro gas turbines, ramp rate constraints are deemed non-binding and are therefore omitted. The operational constraints are solely confined to their power output bounds as follows:

(9)$P_{\mathrm{GT},\min }\le P_{\mathrm{GT},i,t}\le P_{\mathrm{GT},\max }$

Here, $P_{\mathrm{GT},\max }$ and $P_{\mathrm{GT},\min }$ represent the upper and lower bounds of allowable power output for the gas turbine unit, respectively.

3.. Power Output Constraints of PV Units

(10)$0\le P_{\mathrm{PV},t}\le P_{\mathrm{PV},\max }$

Here, $P_{\mathrm{PV},\max }$ denotes the maximum allowable output power of the PV unit.

4.. Power Output Constraints of Wind Turbine Units

(11)$0\le P_{\mathrm{WT},t}\le P_{\mathrm{WT},\max }$

Here, $P_{\mathrm{WT},\max }$ represents the wind turbine unit’s maximum allowable power output, subject to aerodynamic limits and grid compliance requirements.

5.. Operational Constraints of ESSs

The ESS must satisfy power constraints in both charging and discharging directions, which is mathematically defined as follows:

(12)$P_{\mathrm{ESS},i,t}^{\min }\le P_{\mathrm{ESS},i,t}\le P_{\mathrm{ESS},i,t}^{\max }$

where $P_{\mathrm{ESS},i,t}^{\min },P_{\mathrm{ESS},i,t}^{\max }$ respectively represent the upper and lower boundaries of the feasible regulation capacity for the ESS, dynamically allocated to each time slot, t, through day-ahead disaggregation.

6.. Demand Response Operational Constraints

(13)$\left\{\begin{array}{l}0\le P_{\mathrm{IL},i,t}\le {\overline{\delta }}_i^{\mathrm{IL}}{P}_{\mathrm{L},i,t}\\ Q_{\mathrm{IL},i,t}=\left(Q_{\mathrm{L},i,t}/P_{\mathrm{L},i,t}\right)P_{\mathrm{IL},i,t}\end{array}\right.\forall i,\forall t$

Here, ${\overline{\delta }}_i^{\mathrm{IL}}$ denotes the maximum allowable curtailment percentage for interruptible loads, representing the contracted flexibility envelope as a fraction of baseline power.

7.. Linearized Power Flow Constraints

(14)$P_{j,t}={\displaystyle \sum _{i\in a(j)}{P}_{ij,t}}-{\displaystyle \sum _{k\in b(j)}{P}_{jk,t}}$

(15)$Q_{j,t}={\displaystyle \sum _{i\in a(j)}{Q}_{ij,t}}-{\displaystyle \sum _{k\in b(j)}{Q}_{jk,t}}$

(16)$V_{j,t}=V_{i,t}-2(r_{ij}{P}_{ij,t}+x_{ij}{Q}_{ij,t})$

Here, $a(j)$ and $b(j)$ denote the set of terminal nodes of branches, with node $j$ as the starting node, and the set of starting nodes of branches, with node $j$ as the terminal node, respectively. $V_{i,t}$ is the square of the voltage at node i, $V_{j,t}$ is the square of the voltage at node $j$, $P_{ij,t}$ is the active power of line $ij$ at time t, and $Q_{ij,t}$ is the reactive power of line $ij$ at time t. $r_{ij}$ is the active power of line $jk$ at time t. $Q_{jk,t}$ is the reactive power of line $jk$ at time t. $r_{ij}$ is the resistance of line $jk$, and $x_{ij}$ is the reactance of the line ij.

8.. Node Voltage Constraints

(17)$U_{\min }^2<V_{j,t}<U_{\max }^2$

Here, U_min and U_max denote the lower and upper limits of the nodal voltage magnitude, which are respectively set to 0.95 per unit (p.u.) and 1.05 p.u. according to international standards.

9.. Power Line Ampacity Constraints

(18)$P_{ij,t}^2+Q_{ij,t}^2\le S_{ij\max }^2$

Here, S_ijmax denotes the maximum allowable active power flow through line ij, which is constrained by thermal and stability limits.

To enhance computational tractability, the nonlinear constraint (18) is linearized using the method proposed in [25], resulting in the following linear formulation:

(19)$\left\{\begin{array}{l}-S_{ij\max }\le P_{ij,t}\le S_{ij\max }\\ -S_{ij\max }\le Q_{ij,t}\le S_{ij\max }\\ -\sqrt{2}\cdot S_{ij\max }\le P_{ij,t}+Q_{ij,t}\le \sqrt{2}\cdot S_{ij\max }\\ -\sqrt{2}\cdot S_{ij\max }\le P_{ij,t}-Q_{ij,t}\le \sqrt{2}\cdot S_{ij\max }\end{array}\right.$

3.3. Multi-Parametric Optimization Model for VPP Dispatch Instruction Disaggregation

Treating the upper-grid dispatch command as parametric terms in the constraints, we derive a compact multi-parametric optimization model, as follows:

(20)$\begin{array}{l}\underset{x}{\min }F=c^{T}x+d\\ \\ \mathrm{s}.\mathrm{t}. Ax\le b+H\theta \\ \hfill Gx=n+M\theta \end{array}$

where x, c, d, $\mathsf{\theta }$, $\mathbf{A}$, $\mathbf{b}$, $\mathbf{G}$, $\mathbf{H}$, $\mathbf{M}$, and $\mathbf{n}$ all correspond to the variables and coefficients in Equations (1)–(19) above, $\mathbf{x}$ is the decision variable vector, including the power output of micro gas turbines, photovoltaic units, and wind turbines, the charging–discharging power of energy storage systems, and the power of curtailable loads, as well as voltages and branch capacities, $\mathbf{c}$ is the cost coefficient vector, $\mathbf{d}$ is the constant term of the cost, $\mathbf{A}$ is the coefficient matrix of decision variables in all inequality constraints, $\mathsf{\theta }$ is the parameter vector representing the real-time dispatch instructions issued by the superior power grid to the virtual power plant, which means that $P_{\mathrm{Instru} t}$ in (2) is equal to $\theta$, $\mathbf{b}$ is the constant term vector of the right-hand side of all inequality constraints, excluding the dispatch instructions, $\mathbf{H}$ is the coefficient matrix of the parameter vector in all inequality constraints, $\mathbf{n}$ is the constant term vector of the right-hand side of all equality constraints, excluding the dispatch instructions, and $\mathbf{M}$ is the coefficient matrix of the parameter vector in all equality constraints.

4. Computational Framework for Multi-Parametric Program Transformation and Solution

4.1. Critical Region Partitioning and Mapping Relationship Construction

Identification of Active and Inactive Constraints

To determine the initial critical region, $C{R}_0$, the constraint sets must be partitioned into active and inactive constraints. An initial parameter vector, ${\theta }_0$, is specified. Assuming that $x^{*}$ is the optimal solution at ${\theta }_0$, the associated optimization problem is solved to delineate the boundaries of the feasible region.

For any given parameter vector, ${\theta }_0\in M^{*}$, the original multi-parametric programming problem is solved to obtain the optimal solution, $x^{*}$. The sets of active and inactive constraints are subsequently determined as follows:

(21)$A=\left\{i\left|A_i{x}^{*}=b_i+H_i{\theta }_0\right.\right\}$

(22)$A_N\left(R_x\theta +r_x\right)\le b_N+H_N\theta$

2.. Determination of Critical Regions

Following the aforementioned transformations, the expression for the critical region is obtained as follows:

(23)$C{R}_A=\left\{\theta \left|\left(A_N{R}_x-H_N\right)\theta \le b_N-A_N{r}_x,R_{\lambda }\theta +r_{\lambda }\ge 0,\right.\theta \in M^{*}\right\}$

3.. KKT Optimality Conditions

Based on the identified sets of active and inactive constraints, the KKT optimality conditions for the multi-parametric programming problem are formulated below. The Lagrangian function for the primal multi-parametric programming problem is expressed as follows:

(24)$L\left(x,\lambda ,\theta \right)=c^{T}x+d+{\lambda }^T\left(Ax-b-H\theta \right)$

Stationarity condition

(25)$c+A^T\lambda =0$

Complementary slackness condition

(26)$\begin{array}{l}{\lambda }_i\left(A_{i}x-b_i-H_i\theta \right)=0\\ {\lambda }_i\ge 0\end{array}$

4.. Solution of KKT Conditions to Derive Mapping Relationship

The KKT conditions are solved to derive the mapping relationship between the optimal dispatch allocation scheme and the scheduling instructions.

Combining the active constraints with the aforementioned stationarity condition yields the following equations:

(27)$\left[\begin{array}{cc}0& A_A^T\\ A_A^T& 0\end{array}\right]\left[\begin{array}{c}x\\ {\lambda }_A\end{array}\right]=\left[\begin{array}{c}-c\\ b_A+H_A\theta \end{array}\right]$

(28)$\left[\begin{array}{c}x\\ {\lambda }_A\end{array}\right]={\left[\begin{array}{cc}0& A_A^T\\ A_A^T& 0\end{array}\right]}^{-1}\left[\begin{array}{c}-c\\ b_A+H_A\theta \end{array}\right]$

(29)$\left[\begin{array}{cc}{M}_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right]={\left[\begin{array}{cc}0& A_A^T\\ A_A^T& 0\end{array}\right]}^{-1}$

The solution yields the following equations:

(30)$x^{*}\left(\theta \right)=-c{M}_{11}+M_{12}{b}_A+M_{12}{H}_A\theta$

(31)${\lambda }_A^{*}\left(\theta \right)=-c{M}_{21}+M_{22}{b}_A+M_{22}{H}_A\theta$

5.. Derivation of the Mapping Between Optimal Cost and Dispatch Instructions

Under fixed parametric conditions, the multi-parametric programming problem can be transformed into a linear programming problem. Based on the strong duality principle of linear programming, the optimal values of the primal and dual objective functions are equal. This implies that the VPP’s cost function can be represented by the dual problem’s objective function.

Analysis of the derivation process for parametric expressions of critical region boundaries reveals that within each critical region CR, all parameters correspond to identical sets of active and inactive constraints, and the optimal solution to the dual problem is unique within the critical region, CR. Therefore, the dual model is adopted to derive the VPP’s cost function with respect to the parameter vector, $\theta$. The dual problem of the primal multi-parametric programming problem is given as follows:

(32)$\begin{array}{l}\underset{\lambda \ge 0}{\max }={\left(b+H\theta \right)}^T \lambda +d\\ \mathrm{s}.\mathrm{t}. A^T\lambda \ge c \lambda \ge 0\end{array}$

Given the optimal solution, ${\lambda }_0^{*}$, to the dual problem under the parameter vector, ${\theta }_0$, the objective function equation with respect to ${\theta }_0$ within the critical region, $C{R}_0$, is obtained as follows:

(33)$z_0^{*}\left(\theta \right)={\left(b+H\theta \right)}^T {\lambda }_0^{*}+d$

4.2. Recursive Partitioning of Parameter Space

After determining the first critical region, $C{R}_0$ (i.e., identifying the optimal active constraint set under a given parameter and obtaining $C{R}_0$ accordingly), the remaining feasible parameter space must be recursively partitioned so that the entire parameter space is divided into multiple critical regions. The specific procedure is as follows:

Step 1: Generate candidate regions.

Starting from the boundary of the initial critical region, $C{R}_0$, for each constraint that defines the boundary of $C{R}_0$ (in the form $R_i\theta \le r_i$, where this constraint holds as an equality within $C{R}_0$), we take the reverse inequality (i.e., generate $R_i\theta >r_i$) to describe the portion of the parameter space that lies beyond the current critical region boundary. In this way, for each active boundary constraint of $C{R}_0$, we obtain an adjacent “candidate remaining region,” $R_{\mathrm{remain},i}$. These regions represent the parts of the parameter space that extend beyond $C{R}_0$ and have not yet been explored.

Step 2: Determine a new critical region.

For each newly generated candidate region, $R_{\mathrm{remain},i}$, select any point, ${\theta }_{\mathrm{new}}$, within that region (commonly chosen slightly inside the boundary or near the center of the region) as a representative, and solve the original optimization problem once (or equivalently, apply the KKT conditions at that point). Through this single solution, we can identify the optimal active constraint set that corresponds to ${\theta }_{\mathrm{new}}$. If this active set has not appeared in the existing collection of critical regions, then a new critical region, $C{R}_i$, can be defined based on this set of constraints. In other words, $C{R}_i$ is determined by the joint effect of these active constraints, and its range can be obtained by solving the equality and inequality relations formed by these constraints. Thus, we pinpoint a new critical region, $C{R}_i$, within the candidate region.

Step 3: Recursive iteration.

Incorporate the newly discovered critical region, $C{R}_i$, into the set of identified regions, and remove the portion covered by $C{R}_i$ from the remaining space to be partitioned. Next, take $C{R}_i$ as a new starting point and repeat Steps 1 and 2: examine each active boundary constraint of $C{R}_i$, generate further candidate regions, and attempt to discover the next critical region. This process proceeds recursively—either in depth-first or breadth-first fashion—until all uncovered portions of the feasible parameter space have been partitioned into critical regions. During recursion, it is necessary to avoid duplicates; if the active set solved in a candidate region corresponds to a critical region that has already been found, skip that region and continue exploring elsewhere.

Step 4: Complete the partition.

When the entire feasible parameter space has been partitioned into critical regions, the recursive process terminates. At this point, the initial region, $C{R}_0$, and the subsequently found regions, $C{R}_1,C{R}_2,\dots ,C{R}_n$, jointly constitute a partition of the entire parameter space. In other words, the union of these disjoint critical regions covers the feasible parameter domain,$R_{\mathrm{feasible}}$, and the active constraint combination remains unchanged within each region. We denote the union of the remaining partitioned regions as $R_{\mathrm{remain}}=C{R}_1\cup C{R}_2\cup \dots \cup C{R}_n$ (excluding the already partitioned $C{R}_0$) to describe this partitioning result.

Through the recursive partitioning process described above, we eventually obtain an optimal partition of the parameter space. Under this partition, different critical regions correspond to different constraint activation patterns; therefore, the optimal decisions within each region (such as the optimal outputs of resources and the dispatch cost in this problem) can be expressed as affine functions of the original parameters. Consequently, the optimal solution is mapped over the entire parameter space in a piecewise-affine form. Specifically:

For the optimal outputs of VPP resources, the piecewise function can be written as follows:

(34)$x(\theta )=\left\{\begin{array}{ll}{M}^{(0)}\theta +b^{(0)},\hfill & \theta \in C{R}_0,\hfill \\ M^{(1)}\theta +b^{(1)},\hfill & \theta \in C{R}_1,\hfill \\ ⋮\hfill \\ M^{(n)}\theta +b^{(n)},\hfill & \theta \in C{R}_n,\hfill \end{array}\right.$

where the matrix, $M^{(i)}$, and the vector, $b^{(i)}$, are the analytical coefficients for the region $C{R}_i$, derived from the KKT conditions within that region.

Similarly, the optimal cost of the VPP dispatch can be expressed as a piecewise-affine function, as follows:

(35)$z(\theta )=\left\{\begin{array}{ll}{\lambda }^{(0)T}\theta +d^{(0)},\hfill & \theta \in C{R}_0,\hfill \\ {\lambda }^{(1)T}\theta +d^{(1)},\hfill & \theta \in C{R}_1,\hfill \\ ⋮\hfill \\ {\lambda }^{(n)T}\theta +d^{(n)},\hfill & \theta \in C{R}_n,\hfill \end{array}\right.$

where ${\lambda }^{(i)}$ and $d^{(i)}$ are the Lagrange multiplier vector and the constant term corresponding to the region $C{R}_i$.

By means of the recursive partitioning described above, we achieve an optimal mapping across the entire parameter space. Whenever a new dispatch instruction, $\theta$, is issued by the upper-level grid, it is sufficient to determine which critical region it falls into and then directly use the corresponding affine function to compute the optimal resource outputs and dispatch cost. This process eliminates the need for re-solving the optimization problem online, greatly improving response speed and demonstrating the practicality of the proposed method for real-time large-scale scheduling.

4.3. Strategies for Dynamic, Unforeseen Market Conditions or Uncertainties in DERs

From a managerial perspective, this paper discusses the interaction between VPPs and ancillary service markets. The method presented here is based on the assumption that a mutual agreement has been made between the VPP and the market: the VPP is required to submit its hourly adjustable range in advance, and the market issues dispatch instructions within that range. Both parties agree on pricing and payment details beforehand. When dynamic and unforeseen changes in DER availability or market conditions occur, the responsibilities of each party are clearly defined. For example, if changes in DER availability or activation costs occur, the VPP must bear the penalty for failing to meet the previously submitted adjustment range or the loss incurred due to increased DER activation costs. If the dispatch instructions from the ancillary service market exceed the submitted adjustment range, the VPP can refuse to execute.

Building on these managerial aspects, the VPP method proposed in this paper needs to offer technical solutions for handling changes in DER availability or costs. Our method’s “instruction”–“dispatch plan” explicit mapping is designed to be recalibrable when such changes occur. When dynamic, unforeseen changes occur in DER costs, outputs, or available adjustment ranges during the day-ahead stage, it is sufficient to re-trigger the local critical region re-partitioning process for the affected constraints. The previously unaffected regions remain valid. Of course, this still involves an increase in computational effort. We must strike a balance between accuracy and complexity to ensure the robustness of the method. To address this, we adopt a threshold-based approach, with the following strategy:

After executing a dispatch instruction, the latest DER status corresponding to the next dispatch instruction (e.g., 15 min later) is pre-obtained in the ultra-short term. By comparing this with the DER forecast values used in the day-ahead strategy, if the deviation exceeds a certain threshold, ε (e.g., 3% of the day-ahead forecast), we update and adjust the optimization model parameters and re-solve the multi-parametric problem offline to generate the applicable strategy repository for the current period. The 15-min forecast accuracy is typically high enough that further fluctuations in DER resources can be ignored. As long as the strategy repository is updated before the next instruction is issued, the instantaneous response can continue to be guaranteed. The threshold for recalculating the multi-parametric programming model represents the trade-off we need to make between accuracy and computational load.

In addition to the strategy outlined above, in practical applications, considering the 15-min forecasting error and other completely unforeseen parameter changes, we suggest leaving a margin in the dispatch strategy repository (e.g., conservatively estimating key information such as DER available capacity) to enhance robustness against unanticipated fluctuations and ensure execution flexibility. However, this may sacrifice some economic optimality, representing the trade-off between real-time performance and robustness.

5. Case Studies

5.1. Case Study Overview

The case study employs the IEEE 33-node distribution network for simulation analysis, as depicted in Figure 2. The network configuration includes

Two microturbines located at Buses 24 and 6;

One PV unit at Bus 16;

One wind turbine at Bus 32;

Curtailable loads at Bus 20;

ESS at Bus 27.

Key parameters for the microturbines, ESS, and demand response resources are provided in Table 1, Table 2, and Table 3, respectively. Power profiles for PV and wind generation are illustrated in Figure 3. Figure 4 shows the adjustable capacity ranges reported by the VPP. Simulations were conducted in MATLAB R2018a using YALMIP Toolbox with the CPLEX 12.8 solver.

5.2. Critical Region Partitioning and Characterization of Optimal Dispatch Scheme

The mathematical mapping between optimal dispatch costs and scheduling instructions at selected time periods is illustrated in Figure 5. Each color corresponds to one critical region.

Figure 5 demonstrates that at each given time period, the dispatchable space of the system is partitioned into multiple critical regions with well-defined boundaries. Within each critical region, the optimal dispatch cost exhibits a distinct linear characteristic, which is mathematically expressed as a linear function of the scheduling instructions. These linear functions, distributed across critical regions, are combined piecewise to form a complete cost function that characterizes the globally optimal dispatch behavior of the VPP. This piecewise linear representation of the cost function reduces computational complexity and provides an efficient tool for real-time dispatching decisions of the VPP.

Further analysis reveals that at each time period, the resource dispatch scheme and the optimal cost of the VPP are influenced by three factors: the adjustable capacity of DERs, the costs of resources, and network constraints. Taking 12:00 a.m. as an example time period, the critical regions within the adjustable range are partitioned. On this basis, the mapping relationships between dispatch instructions and both (i) resource dispatch schemes and (ii) optimal dispatch costs within each critical region are presented in Table 4.

Table 4 delineates the mapping relationships between scheduling instructions and the resource dispatch schemes/optimal costs across seven distinct critical regions. Within each critical region, the dispatch cost, z, can be expressed as a linear function of the scheduling instructions. At this stage, the partitioning of the dispatch instruction’s feasible region and the visualization of the VPP’s optimal cost function are depicted in Figure 6 and Figure 7, respectively.

At 12:00 a.m., the feasible region for dispatch instructions is defined within the following range: 1671 kW ≤ θ ≤ 6354 kW. Partitioning this parametric feasible region yields seven critical regions. Each critical region corresponds to distinct resource portfolios and operational constraints, reflecting the cost performance attributes of different resource combinations in the VPP. Notably, when the dispatch instruction coincides with the planned generation output, the VPP dispatch cost is non-zero. This is associated with the original dispatch scheme, including PV curtailment.

It can be seen that the optimal dispatch cost increases as the dispatch instruction deviates from the planned generation output. Within the interval 2472 kW ≤ θ ≤ 2639 kW, the cost function exhibits a flatter slope, indicating that the VPP primarily utilizes curtailable loads—the resource with the lowest unit dispatch cost—for regulation. In high-dispatch instruction regions, the cost function slope increases significantly, reflecting the activation of higher-cost reserve resources that elevate marginal costs.

Solving an optimization model with nonlinear network constraints, such as power line ampacity limits, can indeed yield accurate and truly optimal dispatch solutions. However, classical nonlinear optimal power flow (OPF) formulations are typically nonconvex, and solving them requires nonlinear programming techniques with substantial computational complexity and time cost.

For real-time scheduling, especially in the context of a VPP that promptly responds to upper-level grid instructions, the priority is to quickly mobilize DERs so that their aggregate response meets the specified power level. Linearizing the nonconvex model is a common acceleration strategy. Nevertheless, even conventional linearized formulations cannot guarantee instantaneous solution generation upon receiving a dispatch command. This motivates our adoption of a multi-parametric approach, which represents a demand-driven compromise between solution accuracy and computational speed.

Our method, as well as traditional linearized optimization approaches, introduces no error in fulfilling the upper-grid dispatch instruction itself. However, the linearization of network constraints may lead to slight violations (e.g., marginal overloads) in certain branches. To quantify this, we conducted a case study: at the dispatch instruction θ = 4500 kW (point A in Figure 7), the nonlinear model solved by IPOPT required 15.3 s, whereas our “table lookup” method produced a solution in less than 0.3 s, showing a substantial difference in computational efficiency. Comparing the two dispatch schemes, the largest discrepancy in branch power flow occurred at branch 28, with a relative error of 2.7%, which remains acceptable in practical settings.

In practical applications where constraint compliance must be strictly enforced, the linearized constraints can be tightened (i.e., made more conservative) within the offline model to ensure a feasible margin under a linear approximation.

5.3. Sensitivity Analysis of Resource Dispatch Costs

Figure 8 illustrates the VPP’s optimal operational cost function at 12:00 noon under different resource dispatch cost conditions. We increase the cost of ESS to 0.22 CNY/kWh, while reducing wind and PV curtailment costs to 0.10 CNY/kWh and 0.15 CNY/kWh, respectively.

Figure 8 demonstrates that adjusting resource dispatch costs alters the partitioning of critical regions. Each color corresponds to a different critical region. Specifically, after cost parameter adjustments, the dispatchable resources within each critical region exhibit distinct characteristics compared to their pre-adjustment states. This heterogeneity in resource portfolios directly affects the VPP’s optimal dispatch costs, inducing a global downward trend in the cost function curve.

This trend can be illustrated with two examples. First, the significant reduction in wind and PV curtailment costs prioritizes these renewables for downward regulation when dispatch instructions fall below planned generation. Second, escalated ESS costs shift preference to gas turbine resources for upward regulation beyond the planned output. This operational paradigm effectively circumvents high-cost resource activation, achieving enhanced economic efficiency.

The proposed methodology exhibits robust responsiveness to resource cost parameters, efficiently computing corresponding dispatch schemes and optimal costs and demonstrating superior generalization capability across diverse operational scenarios.

5.4. Impact Analysis of Network Constraints on Optimization Results

By comparing scenarios with and without network constraints, critical regions are partitioned, and mapping relationships within each critical region are analyzed. The optimal dispatch cost function curves for both scenarios are illustrated in Figure 9.

Analysis of the graphical results indicates that when distribution network constraints are disregarded, the VPP’s adjustable capacity expands to [1671, 7154] kW. This expansion occurs because PV units and gas turbine units cannot fully utilize their regulation potential under network-constrained conditions due to operational limit violations. Consequently, VPP dispatch models that neglect network constraints may overestimate actual regulation capabilities, potentially leading to challenges in maintaining real-time power balance and compromising grid security.

Network constraints significantly impact VPP dispatch strategy selection. Under ideal network-unconstrained conditions, the VPP would implement economic merit-order dispatch solely based on the marginal costs of individual resources, eliminating the need for resource portfolio optimization. In actual operations, however, network constraints compel the adoption of more sophisticated dispatch strategies, highlighting the significant impact of physical grid constraints on VPP operational decision-making.

5.5. Scalability Analysis and Strategies of the Proposed Method

Number and Complexity of DERs

From a mechanism perspective, since multi-parametric programming is performed in the day-ahead stage, there is sufficient time for solving the problem. In contrast, the “table lookup” solution in the real-time stage has significant scalability. Given that the one-dimensional range is relatively simple and the critical regions are sequentially arranged, as the number and complexity of DERs increase, the proposed method results in an expanded one-dimensional space partitioning. However, the optimization model dimensions and solving difficulty in traditional methods increase exponentially with the number of DERs. Clearly, the “table lookup” solution, which does not involve the “optimization” process, is always simpler than traditional optimization methods. While the proposed method has its inherent scalability limits, these limits are expected to be greater than those of traditional methods. The specific limit is influenced by computational performance.

For example, in the case of the 33-node system described in the paper, with 30 DERs connected to different nodes (multiple DERs at the same node can be simplified and treated as a single DER using aggregation techniques), the number of critical regions obtained is approximately 9000. For real-time scheduling simulation, the lookup time is less than 1 s, with no issues arising.

Of course, the scale of multi-parametric programming in the day-ahead stage significantly increases as the number and complexity of DERs grow, necessitating further analysis and control. The theoretical upper bound on the growth of critical regions is related to the number of active constraint combinations, but in practice, the active sets encountered are much fewer than the combinatorial upper bound (due to the presence of redundant and mutually exclusive constraints). Therefore, the following suppression and compression techniques can be applied:

Redundancy detection: Linear dependence and support set detection are performed on newly generated candidate active sets. If the linear mapping of the new region is consistent with the affine function of an existing region and the boundary set does not expand the feasible domain, the regions can be merged.

ε-merging: In neighboring regions with affine cost function parameters (α, β), a distance < ε can be aggregated to form an approximate piecewise representation, thus controlling the number of critical regions.

Scale of the System

To verify the scalability of the proposed method, the IEEE 123-bus system [26] is adopted. Two ESS units (with a maximum charge/discharge power of 300 kW and an energy capacity of 600 kWh each) were installed at Buses 18 and 67. Two microturbines with a maximum power output of 1000 kW were placed at Buses 25 and 86. Four 500 kW renewable generation (REG) units (two PVs and two WTs) were integrated at Buses 30, 75, 96, and 114. The corresponding results are shown in Figure 10.

Offline stage: The recursive critical region partitioning time increased from 23.4 s to 129.7 s, but this remains easily achievable within the day-ahead window.

Online stage: The table lookup time grew linearly with the number of regions, from 0.21 s to 0.94 s, still meeting the instantaneous response requirement after the instruction is issued.

The number of critical regions increased from 8 to 46.

Figure 10

Optimal cost function of the VPP at 12:00 a.m. in the 123-Bus system.

[Figure omitted. See PDF]

Regarding scalability, the proposed method remains applicable in larger and more complex distribution networks in principle. Although the total number of constraints and the number of possible critical regions will increase with the addition of more nodes and DERs, two factors ensure that the computational demand after scaling can be managed:

As mentioned earlier, each DER or network constraint corresponds to at most one primary region boundary. The number of critical regions is expected to grow approximately linearly with the number of active constraints, which prevents an explosive partitioning of the parameter space.

The parameter dimension in our dispatch disaggregation problem is typically limited (in this study, the upper-level dispatch power is a single parameter). Even with hundreds or thousands of DERs, the required partitioning is of a one-dimensional instruction range, which significantly reduces the complexity of region partitioning compared with high-dimensional cases.

These factors together ensure that the scalability of the method can be efficiently managed, even in larger and more complex distribution network scenarios.

6. Conclusions

This paper proposes a multi-parametric programming-based VPP dispatch instruction disaggregation approach, addressing the computational complexity bottleneck of scheduling methods when massive distributed resources are widely integrated. This achieves the efficient mapping of real-time up-grid dispatch instructions into resource dispatch schemes and optimal VPP scheduling costs. The conclusions of this paper are as follows:

(1) A multi-parametric programming model for VPP dispatch is constructed with up-grid dispatch instructions as parameters. By recursively partitioning the parameter space and solving KKT conditions, expressions for optimal dispatch strategies and cost functions within critical regions are derived, achieving real-time mapping from instruction input to dispatch strategies. This approach generates a parameterized dispatch strategy repository, enabling real-time optimal solutions through simple table lookup during disaggregation and effectively resolving the real-time computational bottleneck in scheduling massive distributed resources.

(2) Comparative case studies analyze the optimal dispatch cost curves of the VPP under varying resource scheduling costs. The results demonstrate that the multi-parametric programming approach dynamically adjusts critical region partitioning and resource dispatch priorities according to heterogeneous operational costs of distributed resources. Therefore, the effectiveness of the proposed instruction disaggregation method for VPPs is verified.

Author Contributions

Conceptualization, Z.Z. and Y.W.; Methodology, Z.Z. and Y.W.; Software, Z.Z.; Validation, Z.Z.; Formal analysis, Y.W.; Data curation, Y.W.; Writing—original draft, Z.Z.; Writing—review & editing, Y.W.; Visualization, Z.Z.; Supervision, Y.W. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Real-time stage virtual power plant instruction disaggregation procedure.

Figure 2 IEEE 33-Bus distribution test system.

Figure 3 Wind and PV power output profiles for a typical day.

Figure 4 Adjustable capacity range curves of dispatchable resources.

Figure 5 Functional mapping between optimal dispatch costs and scheduling instructions at discrete time periods.

Figure 6 Partitioning diagram of the dispatch instruction feasible region at 12:00 a.m.

Figure 7 Optimal cost function of the VPP at 12:00 a.m.

Figure 8 Optimal dispatch cost function of the VPP before and after resource cost adjustment.

Figure 9 Optimal cost functions of the VPP with and without network constraints.

Table 1

Microturbine technical specifications.

No.	Bus	Power Generation Cost Coefficient/CNY (kWh)−1	Capacity/kW
1	24	0.19	1000
2	6	0.14	1500

Table 2

Energy storage system parameters.

Bus	Maximum Charge and Discharge Power/kW	Energy Storage Capacity/kWh	Cost Coefficient/CNY (kWh)−1
27	300	600	0.012

Table 3

Demand response resource parameters.

Demand Response Type	Bus	Available Time	Power/kW	Cost Coefficient/CNY (kWh)−1
Curtailable loads	20	9:00–17:00	200	0.032

Table 4

Dispatch costs and resource dispatch schemes across critical regions at 12:00 a.m.

Demand Response Type	CR1	CR2	CR3	CR4	CR5	CR6	CR7
Objective function z	z = 1000.6 − 0.43θ	z = 481.6 − 0.151θ	z = 572.4 − 0.196θ	z = 117.6 − 0.012θ	z = 0.012θ + 56.4	z = 0.162θ − 339.5	z = 0.43θ − 1984.7
Wind power (kW)	888	888	2472 − 0.785θ	1350	1350	1350	1350
PV power (kW)	1050	2018 − 0.52θ	1200	1200	1200	1200	1200
ESS power (kW)	−78	−78	−78	θ − 2550	θ − 2550	89	89
Curtailable loads (kW)	0	0	0	0	0	0	0
GT 1 power (kW)	0	0	0	0	0	0.429θ − 1133.63	1500
GT 2 power (kW)	0	0	0	0	0	0.571θ − 1505.37	2000

References

1. Shaikh, M.S.; Lin, H.; Zheng, G.; Wang, C.; Lin, Y.; Dong, X. Innovative Hybrid Grey Wolf-Particle Swarm Optimization for Calculating Transmission Line Parameters. Heliyon; 2024; 10, e38555. [DOI: https://dx.doi.org/10.1016/j.heliyon.2024.e38555] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39512928]

2. Shaikh, M.S.; Hua, C.; Jatoi, M.A.; Ansari, M.M.; Qader, A.A. Parameter Estimation of AC Transmission Line Considering Different Bundle Conductors Using Flux Linkage Technique. IEEE Can. J. Electr. Comput. Eng.; 2021; 44, pp. 313-320. [DOI: https://dx.doi.org/10.1109/ICJECE.2021.3069143]

3. Zhou, K.; Peng, N.; Yin, H.; Hu, R. Urban virtual power plant operation optimization with incentive-based demand response. Energy; 2023; 282, 128700. [DOI: https://dx.doi.org/10.1016/j.energy.2023.128700]

4. Lin, C.; Hu, B.; Tai, H.-M.; Shao, C.; Xie, K.; Wang, Y. Performance optimization of VPP in fast frequency control ancillary service provision. Appl. Energy; 2024; 376, 124294. [DOI: https://dx.doi.org/10.1016/j.apenergy.2024.124294]

5. Yin, S.; Ai, Q.; Li, J.; Li, D.; Guo, Q. Trading mode design for a virtual power plant based on main-side consortium blockchains. Appl. Energy; 2022; 325, 119932. [DOI: https://dx.doi.org/10.1016/j.apenergy.2022.119932]

6. Chen, Q.; Lyu, R.; Guo, H.; Su, X. Real-time operation strategy of virtual power plants with optimal power disaggregation among heterogeneous resources. Appl. Energy; 2024; 361, 122876. [DOI: https://dx.doi.org/10.1016/j.apenergy.2024.122876]

7. Liu, Y.; Li, Y.; Wang, Y.; Zhu, J.; Gooi, H.B.; Xin, H. Distributed real-time multi-objective control of a virtual power plant in DC distribution systems. IEEE Trans. Power Deliv.; 2022; 37, pp. 1876-1887. [DOI: https://dx.doi.org/10.1109/TPWRD.2021.3099834]

8. Zhang, X.; Yan, L.; Ullah, Z.; Qazi, H.S.; Hasanien, H.M. A Novel Solution Strategy for Scheduling Optimization of Virtual Power Plant Considering Multiple Participants and Peak Energy Market. Renew. Energy; 2025; 250, 123275. [DOI: https://dx.doi.org/10.1016/j.renene.2025.123275]

9. Gholami, K.; Arif, M.T.; Haque, M.E. IGDT-based participation of hybrid renewable-hydrogen virtual power plants in frequency control ancillary services. Int. J. Hydrogen Energy; 2025; 139, pp. 630-645. [DOI: https://dx.doi.org/10.1016/j.ijhydene.2025.05.280]

10. Shenoy, V.; Serna-Torre, P.; Schoenwald, D.; Hidalgo-Gonzalez, P.; Poveda, J.I. Fast frequency regulation of virtual power plants via Droop Reset Integral Control (DRIC). Electr. Power Syst. Res.; 2024; 235, 110762. [DOI: https://dx.doi.org/10.1016/j.epsr.2024.110762]

11. Rahimiyan, M.; Baringo, L. Real-time energy management of a smart virtual power plant. IET Gener. Transm. Distrib.; 2019; 13, pp. 2015-2023. [DOI: https://dx.doi.org/10.1049/iet-gtd.2018.5637]

12. Abdolrasol, M.G.M.; Hannan, M.A.; Mohamed, A.; Amiruldin, U.A.U.; Abidin, I.B.Z.; Uddin, M.N. An optimal scheduling controller for virtual power plant and microgrid integration using the binary backtracking search algorithm. IEEE Trans. Ind. Appl.; 2018; 54, pp. 2834-2844. [DOI: https://dx.doi.org/10.1109/TIA.2018.2797121]

13. Gu, N.; Cui, J.; Wu, C. An auto-tuned robust dispatch strategy for virtual power plants to provide multi-stage real-time balancing service. IEEE Trans. Smart Grid; 2023; 14, pp. 4494-4507. [DOI: https://dx.doi.org/10.1109/TSG.2023.3265398]

14. Mei, Z.; Jiang, X.; Wang, K. Multi-objective coordinated optimal scheduling of virtual power plants based on demand side response. Proceedings of the 2022 IEEE Asia-Pacific Conference on Image Processing, Electronics and Computers (IPEC); Dalian, China, 14–16 April 2022; pp. 753-757.

15. Xu, B.; Luan, W.; Yang, J.; Zhao, B.; Long, C.; Ai, Q.; Xiang, J. Integrated three-stage decentralized scheduling for virtual power plants: A model-assisted multi-agent reinforcement learning method. Appl. Energy; 2024; 376, 123985. [DOI: https://dx.doi.org/10.1016/j.apenergy.2024.123985]

16. Yi, Z.; Xu, Y.; Gu, W.; Wu, W. A multi-time-scale economic scheduling strategy for virtual power plant based on deferrable loads aggregation and disaggregation. IEEE Trans. Sustain. Energy; 2019; 11, pp. 1332-1346. [DOI: https://dx.doi.org/10.1109/TSTE.2019.2924936]

17. Chang, W.; Dong, W.; Wang, Y.; Yang, Q. Two-stage coordinated operation framework for virtual power plant with aggregated multi-stakeholder microgrids in a deregulated electricity market. Renew. Energy; 2022; 199, pp. 943-956. [DOI: https://dx.doi.org/10.1016/j.renene.2022.09.037]

18. Liu, X.; Lin, X.; Qiu, H.; Li, Y.; Huang, T. Optimal aggregation and disaggregation for coordinated operation of virtual power plant with distribution network operator. Appl. Energy; 2024; 376, 124142. [DOI: https://dx.doi.org/10.1016/j.apenergy.2024.124142]

19. Huang, J.; Li, H.; Zhang, Z. Review of Virtual Power Plant Response Capability Assessment and Optimization Dispatch. Technologies; 2025; 13, 216. [DOI: https://dx.doi.org/10.3390/technologies13060216]

20. Ahmadi-Moshkenani, P.; Johansen, T.A.; Olaru, S. Combinatorial approach toward multiparametric quadratic programming based on characterizing adjacent critical regions. IEEE Trans. Autom. Control.; 2018; 63, pp. 3221-3231. [DOI: https://dx.doi.org/10.1109/TAC.2018.2791479]

21. Akbari, A.; Barton, P.I. An improved multi-parametric programming algorithm for flux balance analysis of metabolic networks. J. Optim. Theory Appl.; 2018; 178, pp. 502-537. [DOI: https://dx.doi.org/10.1007/s10957-018-1281-x]

22. Burnak, B.; Katz, J.; Pistikopoulos, E.N. A space exploration algorithm for multiparametric programming via Delaunay triangulation. Optim. Eng.; 2020; 22, pp. 555-579. [DOI: https://dx.doi.org/10.1007/s11081-020-09535-6]

23. Arce, A.; Panos, C.; Bordons, C.; Pistikopoulos, E.N. Design and experimental validation of an explicit MPC controller for regulating temperature in PEM fuel cell systems. IFAC Proc. Vol.; 2011; 44, pp. 2476-2481. [DOI: https://dx.doi.org/10.3182/20110828-6-IT-1002.02819]

24. Liu, J.; Wunsch, D.C.; Wang, S.; Bo, R. Multi-parametric analysis for mixed integer linear programming: An application to transmission upgrade and congestion management. Sustain. Energy Grids Netw.; 2024; 40, 101563. [DOI: https://dx.doi.org/10.1016/j.segan.2024.101563]

25. Chen, X.; Wu, W.; Zhang, B.; Lin, C. Data-Driven DG Capacity Assessment Method for Active Distribution Networks. IEEE Trans. Power Syst.; 2017; 32, pp. 3946-3957. [DOI: https://dx.doi.org/10.1109/TPWRS.2016.2633299]

26. Zhang, S.; Ge, S.; Liu, H.; Zhao, B.; Ni, C.; Hou, G.; Wang, C. Region-based flexibility quantification in distribution systems: An analytical approach considering spatio-temporal coupling. Appl. Energy; 2024; 355, 122175. [DOI: https://dx.doi.org/10.1016/j.apenergy.2023.122175]