Dynamic Boundary Dissemination to Virtual Power Plants for Congestion and Voltage Management in Power Distribution Networks

Abstract

Virtual power plants (VPPs) are optimized to maximize profits by efficiently scheduling their resources. However, dynamic power trading over existing distribution networks can lead to voltage disturbances and branch congestion, posing risks to network security. Moreover, distribution network service providers (DNSPs) face the added challenge of managing VPP operations while complying with privacy preservation. To address these challenges, this paper proposes a decentralized co-optimization technique for integrating VPPs into distribution networks. The approach enables DNSPs to define dynamic operational boundaries for VPPs, effectively mitigating network congestion and voltage fluctuations while ensuring privacy. Additionally, the proposed convex optimization framework allows the publication of operational boundaries for multiple VPPs with minimal computational effort, making it suitable for real-time applications. The effectiveness of the technique is validated using the IEEE benchmark network connected with electricity–hydrogen VPPs. Results demonstrate that the proposed approach maintains voltage levels within standard limits and prevents branch congestion, confirming its suitability for stable and secure grid operations.

1. Introduction

The widespread deployment of distributed energy resources (DERs), such as photovoltaic (PV) systems and energy storage, presents valuable opportunities to address environmental concerns and price volatility [1,2]. However, the high penetration of these resources introduces significant challenges for distribution networks, which have traditionally been designed for unidirectional power flow from the power plants to the customers [3,4]. The bidirectional flow introduced by DERs can result in voltage instability and congestion, potentially exceeding network capacity. Previous studies [5,6] explored direct management of behind-the-meter DERs by DNSPs. However, this approach is impractical in developed countries like Australia due to several factors. Firstly, the vast number of DERs makes individual management by DNSPs infeasible. Secondly, privacy regulations restrict DNSPs from accessing DER data across various power system layers [7].

To overcome these challenges, the concept of VPPs [8,9] has been introduced to aggregate DERs, allowing a local agent to manage them and interact with the broader power system as a large-scale entity capable of buying or selling power. While forming VPPs offers significant technical and financial benefits, it also introduces challenges, particularly at the distribution level, such as exceeding branch capacity limitations. Since distribution networks were originally designed for unidirectional power flow, uncoordinated power trading with VPPs can cause congestion and voltage instability and may exceed branch capacity limits. Without proper coordination between VPPs and DNSPs, these challenges could lead to severe disruptions, including potential system-wide blackouts [10]. To mitigate these risks, various strategies have been proposed in the literature and are discussed in the following sections.

Several methods have been proposed for integrating VPPs and managing congestion in distribution networks. For instance, a congestion relief framework introduced in [11] enables DNSPs to interact with VPP service providers through a two-step iterative process. Initially, DNSPs assess VPP bids and perform power flow analysis to detect congestion. If congestion is identified, the network is reconfigured to resolve the issue. If reconfiguration proves insufficient, VPPs are signaled to reschedule their resources. This process repeats until congestion is eliminated. While effective, the approach has limitations, including the use of a nonlinear formulation and an iterative Newton–Raphson method, which imposes a high computational burden, making it impractical for real-time electricity markets. Additionally, network reconfiguration can introduce power quality issues [12]. Another strategy proposed in [13] employs a two-stage congestion relief method using cascade hydro–PV–pumped storage VPPs. The first stage minimizes congestion costs, followed by VPP power corrections to optimize three objectives using genetic algorithms (GAs) and an enhanced non-dominated GA (NSGA-II). The results demonstrated cost reductions through VPP cooperation, with the enhanced NSGA-II outperforming previous versions. However, the nonlinear formulation remains computationally expensive, limiting its feasibility in real-time market operations.

To address these challenges, convex optimization has been identified as a preferable alternative due to its guaranteed convergence within a finite time, making it suitable for dynamic scenarios where DNSPs must promptly publish operational boundaries while ensuring network stability. Recognizing this need, the researchers in [14] proposed a convex bi-level optimization model for VPP integration. Although this approach significantly improved congestion management, it raised privacy concerns, as the use of Karush–Kuhn–Tucker (KKT) conditions to simplify the bi-level problem into a single-level formulation potentially exposes sensitive service provider information.

While the aforementioned studies have made significant contributions to managing voltage and congestion in distribution networks with VPP integration, most face challenges related to data sharing among service providers, which may not be practical in evolving energy markets. To address this issue and ensure privacy while maintaining voltage stability, researchers have explored pricing mechanisms such as Distribution Locational Marginal Pricing (DLMP) [15,16,17] and transactive energy models [18,19]. These methods offer the advantage of clearly delineating the roles of DNSPs. However, their practical implementation—such as enabling DNSPs to issue price signals—faces regulatory challenges due to the transition from traditional passive distribution networks to modern active networks. In most cases, DNSPs do not have the authority to issue price signals, as this responsibility typically falls to independent market operators, such as the Australian Energy Market Operator (AEMO) [7]. In particular, the DNSPs are facing significant challenges in congestion and voltage fluctuations with the increase of bi-directional trading and high penetration of DERs [20]. Traditional approaches to these issues rely on centralized or iterative techniques that may not be able to cope with the dynamic nature of modern energy markets. Consequently, the distribution networks face branch congestion [21] and voltage fluctuations [3] driven by the dynamic power trading of VPPs with markets [22]. For instance, each branch of a distribution network has a designed capacity for power transfer, usually determined by its cable capacity. If the power trading is not dynamically managed, then secure network operation becomes infeasible [23].

These challenges highlight the need for a coordinated strategy that allows DNSPs to manage power trading with VPPs while maintaining technical standards and ensuring secure operations. However, DNSPs typically lack direct access to VPP data [24], necessitating the development of a framework that guarantees data privacy. The key limitations of existing methods are summarized as follows:

• Confidentiality concerns: Conventional methods compromise the privacy of both DNSPs and VPPs due to centralized data sharing.
• Lack of integrated solutions: Current techniques fail to jointly address congestion and voltage control in distribution networks.
• Slow response times: Traditional approaches cannot dynamically publish operational boundaries for VPPs quickly enough.
• Computational inefficiency: Dependency on iterative algorithms makes them unsuitable for real-time applications

To overcome these challenges, this work proposes a novel co-optimization framework that facilitates dynamic collaboration between DNSPs and VPPs. The framework ensures secure network operation while preserving the privacy of all stakeholders. In essence, DNSPs dynamically publish available capacities, allowing VPPs to operate within safe boundaries, thereby mitigating congestion and voltage instability risks. The unique contributions of this approach are as follows:

• Dynamic operational boundaries: Unlike conventional congestion mitigation and voltage control methods, the proposed co-optimization framework dynamically allocates operational limits to VPPs, effectively addressing congestion and voltage fluctuations.
• Decentralized data privacy: The framework ensures minimal data exchange, sharing only permissible operational boundaries set by DNSPs and bids from VPPs, aligning with emerging regulatory requirements for unbundled electricity grids.
• Technology integration: Advanced features such as power-to-hydrogen systems, vehicle-to-grid integration, and demand response mechanisms are incorporated, ensuring the framework’s adaptability to diverse DER technologies.
• Efficient optimization: A convex optimization model is formulated to achieve globally optimal solutions efficiently using commercial solvers like GUROBI. This framework provides a robust, scalable, and privacy-preserving solution for integrating VPPs into modern distribution networks

The remainder of this paper is organized as follows: Section 2 presents the modeling of VPPs, while Section 3 details the proposed approach for publishing operational boundaries. Section 4 discusses the results, followed by the conclusion in Section 5.

2. VPP Component Modeling

This section elaborates on the component modeling of VPPs.

2.1. Objective Function

The primary objective of the VPP is to maximize its profit, which is defined as the revenue generated from power sales to the market minus both internal operational expenditures (IOPEX) and external operational expenditures (EOPEX), as shown in Equation (1). Revenue, as expressed in Equation (2), is derived from multiple sources, including market power sales, electric vehicle aggregation (EVA) charging, hydrogen selling to hydrogen fuel cell vehicles (HFCVs), and power sold to customers. Equation (3) has 6 components, including the operation costs of dispatchable generators (DGs), microturbines (MTs), battery storage systems (BSSs), and the maintenance costs of electrolyzers, photovoltaic systems, and wind turbines. The variable costs of other devices, such as electrolyzers, are considered in the power balance given by Equations (4) and (5). The expression for EOPEX, as explained in Equation (4), includes the cost of buying power from the market and payments to EVA owners, represented by the first and second terms, respectively.

$$\mathrm{O}\mathrm{F}=\max \left[\mathrm{R}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{u}\mathrm{e}-\mathrm{I}\mathrm{O}\mathrm{P}\mathrm{E}\mathrm{X}-\mathrm{E}\mathrm{O}\mathrm{P}\mathrm{E}\mathrm{X}\right]$$

(1)

$$\begin{gathered} \mathrm{R}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{u}\mathrm{e}=\sum _{\mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}}{\mathsf{\pi }}^{\mathrm{s}}\left(\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{{\mathrm{C}}_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{s},\mathrm{t}}\mathrm{P}}_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{s},\mathrm{t}}+\sum _{\mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{s}\mathrm{e}\mathrm{l}\mathrm{l}}^{\mathrm{s},\mathrm{t}}{\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ +\sum _{\mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{f}\mathrm{c}\mathrm{v}}^{\mathrm{h},\mathrm{s}}{\mathrm{H}}_{\mathrm{f}\mathrm{c}\mathrm{v}}^{\mathrm{h},\mathrm{s},\mathrm{t}}+\sum _{\mathrm{c}\in {\Omega }_{\mathrm{c}\mathrm{m}\mathrm{r}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{c}\mathrm{m}\mathrm{r}}^{\mathrm{c}}{\mathrm{P}}_{\mathrm{c}\mathrm{m}\mathrm{r}}^{\mathrm{c},\mathrm{s},\mathrm{t}}\right) \end{gathered}$$

(2)

$$\begin{gathered} \mathrm{I}\mathrm{O}\mathrm{P}\mathrm{E}\mathrm{X}=\sum _{\mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}}{\mathsf{\pi }}^{\mathrm{s}}\left(\sum _{\mathrm{d}\in {\Omega }_{\mathrm{d}\mathrm{g}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d}}{\mathrm{P}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{s},\mathrm{t}} +\sum _{\mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{O}\&\mathrm{M},\mathrm{M}\mathrm{T}}^{\mathrm{h}}{\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}} \\ \hspace{1em}+\sum _{\mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{c}\mathrm{y}\mathrm{c},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b}}\left|{\mathrm{B}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}-{\mathrm{B}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}-1}\right| \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\sum _{\mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{O}\&\mathrm{M},\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h}}{\mathrm{H}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{s},\mathrm{t}}+\sum _{\mathrm{p}\in {\Omega }_{\mathrm{P}\mathrm{V}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{O}\&\mathrm{M},\mathrm{P}\mathrm{V}}^{\mathrm{p}}{\mathrm{P}}_{\mathrm{P}\mathrm{V}}^{\mathrm{p},\mathrm{s},\mathrm{t}} \\ +\sum _{\mathsf{\omega }\in {\Omega }_{\mathrm{W}\mathrm{T}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{O}\&\mathrm{M},\mathrm{W}\mathrm{T}}^{\mathsf{\omega }}{\mathrm{P}}_{\mathrm{W}\mathrm{T}}^{\mathsf{\omega },\mathrm{s},\mathrm{t}}\right) \end{gathered}$$

(3)

$$\mathrm{E}\mathrm{O}\mathrm{P}\mathrm{E}\mathrm{X}=\sum _{\mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}}{\mathsf{\pi }}^{\mathrm{s}}\left(\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{{\mathrm{C}}_{\mathrm{b}\mathrm{u}\mathrm{y}}^{\mathrm{s},\mathrm{t}}\mathrm{P}}_{\mathrm{b}\mathrm{u}\mathrm{y}}^{\mathrm{s},\mathrm{t}}+\sum _{\mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}\left({{\mathrm{C}}_{\mathrm{r}\mathrm{e}\mathrm{m}}^{\mathrm{e}}+\mathrm{C}}_{\mathrm{b}\mathrm{u}\mathrm{y}}^{\mathrm{s},\mathrm{t}}\right){\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}\right)$$

(4)

In Equation (3), the absolute term related to battery cycling is nonlinear and requires linearization. To achieve this, a binary auxiliary variable, denoted as ${\mathrm{B}\mathrm{C}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}$, is introduced. This variable replaces the absolute term. Consequently, the linearized formulation of the VPP’s bidding strategy is revised, incorporating the new variables into the objective function, which is subject to the associated constraints as follows.

$$\begin{gathered} \mathrm{I}\mathrm{O}\mathrm{P}\mathrm{E}\mathrm{X}=\sum _{\mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}}{\mathsf{\pi }}^{\mathrm{s}}\left(\sum _{\mathrm{d}\in {\Omega }_{\mathrm{d}\mathrm{g}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d}}{\mathrm{P}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{s},\mathrm{t}} +\sum _{\mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{O}\&\mathrm{M},\mathrm{M}\mathrm{T}}^{\mathrm{h}}{\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\sum _{\mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{c}\mathrm{y}\mathrm{c},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b}}{\mathrm{B}\mathrm{C}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}+\sum _{\mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{O}\&\mathrm{M},\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h}}{\mathrm{H}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{s},\mathrm{t}} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\sum _{\mathrm{p}\in {\Omega }_{\mathrm{P}\mathrm{V}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{O}\&\mathrm{M},\mathrm{P}\mathrm{V}}^{\mathrm{p}}{\mathrm{P}}_{\mathrm{P}\mathrm{V}}^{\mathrm{p},\mathrm{s},\mathrm{t}}+\sum _{\mathsf{\omega }\in {\Omega }_{\mathrm{W}\mathrm{T}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{C}}_{\mathrm{O}\&\mathrm{M},\mathrm{W}\mathrm{T}}^{\mathsf{\omega }}{\mathrm{P}}_{\mathrm{W}\mathrm{T}}^{\mathsf{\omega },\mathrm{s},\mathrm{t}}\right) \end{gathered}$$

(5)

It is notable that this substitution by itself is not enough and the newly introduced auxiliary variables are also constrained as follows:

$${\mathrm{B}\mathrm{C}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}\ge {\mathrm{B}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}-{\mathrm{B}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}-1}, \forall \mathrm{t}>1\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(6)

$${\mathrm{B}\mathrm{C}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}\ge {{\mathrm{B}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}-1}-\mathrm{B}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}>1\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(7)

2.2. Operational Constraints

Here, the constraints related to the VPPs are provided as follows:

2.2.1. Battery System

The state of charge (SoC) of the battery storage systems (BSSs) at time t = 1 is calculated using Equation (8). For subsequent time periods where t > 1, the SoC is determined using Constraint (9). To ensure the SoC of the BSSs is maintained between its capacity during both charging and discharging, Constraint (10) is deployed. Additionally, the charging and discharging processes are constrained by their respective limits, as outlined in Equations (11) and (12). Furthermore, the initial and final SoC of the BSSs are equal via constraints specified in Equation (13).

$$\begin{gathered} {\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}={\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}+{\mathsf{\eta }}_{\mathrm{c}\mathrm{h}, \mathrm{b}\mathrm{a}\mathrm{t}}{\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}}-\left({\displaystyle \frac{1}{{\mathsf{\eta }}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{b}\mathrm{a}\mathrm{t}}}}\right){\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}} , \forall \mathrm{t}=1, \\ \forall \mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}} , \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}} \end{gathered}$$

(8)

$$\begin{gathered} {\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}={\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}-1}+{\mathsf{\eta }}_{\mathrm{c}\mathrm{h},\mathrm{b}\mathrm{a}\mathrm{t}}{\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}}-\left({\displaystyle \frac{1}{{\mathsf{\eta }}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{b}\mathrm{a}\mathrm{t}}}}\right){\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}} , \forall \mathrm{t}>1\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \\ \forall \mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}} \end{gathered}$$

(9)

$${\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}\le {\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(10)

$${\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}\le {\mathrm{B}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}{\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(11)

$${\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}\le \left(1-{\mathrm{B}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}\right){\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(12)

$${\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},{\mathrm{t}}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}}={\mathrm{S}}_{\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},{\mathrm{t}}^{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}}, \mathrm{e}.\mathrm{g}., {\mathrm{t}}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}=0, {\mathrm{t}}^{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}=24,\forall \mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(13)

2.2.2. EV System

In the future, the use of electric vehicles, as a cleaner alternative to fossil fuel-based vehicles will increase [25]. The following constraints enable VPPs to effectively incorporate electric vehicle aggregation (EVA) into their bidding strategy: Constraints (14) and (15) specify the state of charge for times t = 1 and t > 1. Constraint (16) ensures the state of charge remains within their limits, while Constraints (17) and (18) regulate charging and discharging within permissible ranges. Constraint (19) guarantees that the state of charge for ${\mathrm{e}}^{\mathrm{t}\mathrm{h}}$ EVA is restored to its initial value at the final time, enabling it for next day scheduling.

$$\begin{gathered} {\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}={\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}+{\mathsf{\eta }}_{\mathrm{c}\mathrm{h}, \mathrm{e}\mathrm{v}}{\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}}-\left({\displaystyle \frac{1}{{\mathsf{\eta }}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}}}\right){\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}}-\left({\displaystyle \frac{1}{{\mathsf{\eta }}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}}}\right){\mathrm{P}}_{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}} , \forall \mathrm{t}=1, \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \forall \mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}} \end{gathered}$$

(14)

$$\begin{gathered} {\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}={\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}-1}+{\mathsf{\eta }}_{\mathrm{c}\mathrm{h},\mathrm{e}\mathrm{v}}{\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}}-\left({\displaystyle \frac{1}{{\mathsf{\eta }}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}}}\right){\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}} -\left({\displaystyle \frac{1}{{\mathsf{\eta }}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}}}\right){\mathrm{P}}_{\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{p},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}}, \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \forall \mathrm{t}>1\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}} , \forall \mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}} \end{gathered}$$

(15)

$${\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}\le {\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(16)

$${\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}\le {{\mathrm{B}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(17)

$${\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}\le {\left(1-{\mathrm{B}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}\right)\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(18)

$${\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},{\mathrm{t}}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}}={\mathrm{S}}_{\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},{\mathrm{t}}^{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}}, \mathrm{e}.\mathrm{g}., {\mathrm{t}}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}=0, {\mathrm{t}}^{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}=24,\forall \mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(19)

2.2.3. Power-to-Hydrogen System

The concepts of power-to-hydrogen (P2H) and hydrogen-to-power (H2P) in the context of hydrogen refueling stations are mathematically outlined as follows. The state of hydrogen in the of ${\mathrm{h}}^{\mathrm{t}\mathrm{h}}$ storage at t = 1 and t > 1 is respectively handled using Constraints (20) and (21). The P2H concept, represented by the charging of hydrogen storage via electrolyzers, is articulated in Constraint (22). Constraint (23) details the hydrogen discharged to meet the demands of HFCVs and the consumption by microturbines, representing the H2P concept. The electrolyzers’ electricity consumption for hydrogen production is captured in Constraint (24), while Constraint (25) quantifies the hydrogen consumed by microturbines. Equation (26) is the boundary for power generated by MTs. Ramp rate limitations for microturbines are managed by Constraint (27). The state of hydrogen is bound by Constraint (28) and hydrogen discharge boundaries are defined in Equation (29). Constraints (30) and (31) ensure that the hydrogen produced and consumed by the electrolyzers and MTs, respectively, remain within specified limits. Furthermore, Constraint (32) ensures that the initial hydrogen storage level matches the final level at the end of operations. To ensure reliable refueling for HFCVs, a reliability constraint is introduced in Equation (33).

$${\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}={\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}+{\mathrm{H}}_{\mathrm{c}\mathrm{h},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}}-{\mathrm{H}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}} , \forall \mathrm{t}=1, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(20)

$${\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}={\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}-1}+{\mathrm{H}}_{\mathrm{c}\mathrm{h},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}}-{\mathrm{H}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}{\Delta }^{\mathrm{t}} , \forall \mathrm{t}>1\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}},\forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(21)

$${\mathrm{H}}_{\mathrm{c}\mathrm{h},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}={\mathrm{H}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(22)

$${\mathrm{H}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}= {\mathrm{H}}_{\mathrm{F}\mathrm{C}\mathrm{V}}^{\mathrm{h},\mathrm{s},\mathrm{t}}+{\mathrm{H}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(23)

$${\mathrm{P}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{s},\mathrm{t}}= {\displaystyle \frac{{\mathsf{\Psi }}_{\mathrm{e}\mathrm{i}\mathrm{h}}\times {\mathrm{H}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{s},\mathrm{t}}}{{\mathsf{\eta }}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(24)

$${\mathrm{H}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}}={\displaystyle \frac{ {\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}}}{{\mathsf{\eta }}_{\mathrm{M}\mathrm{T}}\times {\mathsf{\Psi }}_{\mathrm{e}\mathrm{i}\mathrm{h}}}} , \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(25)

$${\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}}\le {\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(26)

$${\mathrm{R}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{D}\mathrm{o}\mathrm{w}\mathrm{n}}\le {\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}}-{\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}-1}\le {\mathrm{R}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{U}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(27)

$${\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}\le {\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(28)

$${\mathrm{H}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{H}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}\le {\mathrm{H}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}},\forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(29)

$${\mathrm{H}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{H}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{s},\mathrm{t}}\le {\mathrm{H}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}},\forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(30)

$${\mathrm{H}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{H}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}}\le {\mathrm{H}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}},\forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(31)

$${\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},{\mathrm{t}}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}}={\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},{\mathrm{t}}^{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}}, \mathrm{e}.\mathrm{g} {\mathrm{t}}^{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}}=0, {\mathrm{t}}^{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}=24,\forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(32)

$${\mathrm{S}}_{\mathrm{h}\mathrm{g}\mathrm{n}}^{\mathrm{h},\mathrm{s},\mathrm{t}}\ge \sum _{\widehat{\mathrm{t}}\in {\Omega }_{\mathrm{S}\mathrm{R}\mathrm{e}\mathrm{l}}}{\mathrm{H}}_{\mathrm{F}\mathrm{C}\mathrm{V}}^{\mathrm{h},\mathrm{s},\mathrm{t}+\widehat{\mathrm{t}}}, \forall \mathrm{t}\in \left\{1,\dots ,{\mathrm{t}}^{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}-{\widehat{\mathrm{t}}}^{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}\right\}, \forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(33)

2.2.4. Demand Response

Demand response programs (DRPs) adjust flexible load demands to enhance system efficiency [26]. Constraints (34) and (35) set the limits for shifting active loads, while Constraints (36) and (37) establish the boundaries for reactive loads. Additionally, Constraints (38) and (39) ensure that the reductions in active and reactive load demands are balanced by corresponding increases.

$$0\le {\mathrm{P}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{U}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}^{\mathrm{l},\mathrm{s},\mathrm{t}}\le {\mathrm{P}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}}^{\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{l}\in {\Omega }_{\mathrm{d}\mathrm{r}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(34)

$$0\le {\mathrm{P}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{D}\mathrm{o}\mathrm{w}\mathrm{n}}^{\mathrm{l},\mathrm{s},\mathrm{t}}\le {\mathrm{P}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}}^{\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{l}\in {\Omega }_{\mathrm{d}\mathrm{r}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(35)

$$0\le {\mathrm{Q}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{U}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}^{\mathrm{l},\mathrm{s},\mathrm{t}}\le {\mathrm{Q}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}}^{\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{l}\in {\Omega }_{\mathrm{d}\mathrm{r}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(36)

$$0\le {\mathrm{Q}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{D}\mathrm{o}\mathrm{w}\mathrm{n}}^{\mathrm{l},\mathrm{s},\mathrm{t}}\le {\mathrm{Q}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}}^{\mathrm{l},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{l}\in {\Omega }_{\mathrm{d}\mathrm{r}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(37)

$$\sum _{\mathrm{l}\in {\Omega }_{\mathrm{D}\mathrm{R}\mathrm{P}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{P}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{D}\mathrm{o}\mathrm{w}\mathrm{n}}^{\mathrm{l},\mathrm{s},\mathrm{t}}=\sum _{\mathrm{l}\in {\Omega }_{\mathrm{D}\mathrm{R}\mathrm{P}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{P}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{U}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}^{\mathrm{l},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{l}\in {\Omega }_{\mathrm{d}\mathrm{r}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(38)

$$\sum _{\mathrm{l}\in {\Omega }_{\mathrm{D}\mathrm{R}\mathrm{P}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{Q}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{D}\mathrm{o}\mathrm{w}\mathrm{n}}^{\mathrm{l},\mathrm{s},\mathrm{t}}=\sum _{\mathrm{l}\in {\Omega }_{\mathrm{D}\mathrm{R}\mathrm{P}}}\sum _{\mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{Q}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{U}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}^{\mathrm{l},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{l}\in {\Omega }_{\mathrm{d}\mathrm{r}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(39)

2.2.5. Dispatchable Generator

Dispatchable generators (DGs) must operate within their specified generation capacity, as per Constraint (40). Additionally, Constraint (41) imposes a ramp rate limitation to ensure stable operation and prevent sudden fluctuations in generation.

$${\mathrm{P}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathrm{P}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{s},\mathrm{t}}\le {\mathrm{P}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{m}\mathrm{a}\mathrm{x}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{d}\in {\Omega }_{\mathrm{d}\mathrm{g}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(40)

$${\mathrm{R}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{D}\mathrm{o}\mathrm{w}\mathrm{n}}\le {\mathrm{P}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{s},\mathrm{t}}-{\mathrm{P}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{s},\mathrm{t}-1}\le {\mathrm{R}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{U}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{d}\in {\Omega }_{\mathrm{d}\mathrm{g}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(41)

2.2.6. Power Flow

In this subsection, we present the constraints related to the internal network of VPPs to incorporate the technical aspects of VPPs into consideration. Constraints (42) and (43) are utilized to determine the active and reactive power flows in the VPP local network, respectively. The voltage of buses is calculated using Constraint (44). Constraints (45) and (46) respectively ensure the balance of active and reactive powers within the VPP’s network. Moreover, Constraint (47) maintains voltage magnitudes within acceptable levels, while Constraint (48) ensures that the power flowing through each branch does not surpass its maximum boundary.

$${\mathrm{P}}_{\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}=\sum _{\mathrm{u}\in {\Omega }_{{\mathrm{p}\mathrm{a}\mathrm{r}}_{\mathrm{v}}}}{\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v},\mathrm{s},\mathrm{t}}-\sum _{\mathrm{w}\in {\Omega }_{{\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{l}}_{\mathrm{v}}}}{\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v}\mathrm{w},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{v}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{V}\mathrm{p}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(42)

$${\mathrm{Q}}_{\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}=\sum _{\mathrm{u}\in {\Omega }_{{\mathrm{p}\mathrm{a}\mathrm{r}}_{\mathrm{v}}}}{\mathrm{Q}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v},\mathrm{s},\mathrm{t}}-\sum _{\mathrm{w}\in {\Omega }_{{\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{l}}_{\mathrm{v}}}}{\mathrm{Q}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v}\mathrm{w},\mathrm{s},\mathrm{t}},\forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{v}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{V}\mathrm{p}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(43)

$$\begin{gathered} {\mathrm{V}}_{\mathrm{b},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}\approx {\mathrm{V}}_{\mathrm{b},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u},\mathrm{s},\mathrm{t}}-{\displaystyle \frac{\left({\mathrm{R}}_{\mathrm{L},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v}}{\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v},\mathrm{s},\mathrm{t}}+{\mathrm{X}}_{\mathrm{L},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v}}{\mathrm{Q}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v},\mathrm{s},\mathrm{t}}\right)}{{\mathrm{V}}_{{\mathrm{s}}_0,\mathrm{V}\mathrm{p}\mathrm{p}}}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{u},\mathrm{v} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{V}\mathrm{P}\mathrm{P}}, {\mathrm{s}}_0 \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{V}\mathrm{P}\mathrm{P} \mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t} \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} , \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}} \end{gathered}$$

(44)

For bus v which is connected with DG d, PV p, WT $\mathsf{\omega }$, BSS b, EVA e, hydrogen system $\mathrm{h}$, and DRP $\mathrm{l}$, $\mathrm{d}\in {\Omega }_{\mathrm{d}\mathrm{g}},\mathrm{p}\in {\Omega }_{\mathrm{P}\mathrm{V}},\mathsf{\omega }\in {\Omega }_{\mathrm{W}\mathrm{T}},\mathrm{b}\in {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{t}},\mathrm{e}\in {\Omega }_{\mathrm{e}\mathrm{v}}$, $\forall \mathrm{h}\in {\Omega }_{\mathrm{h}\mathrm{g}\mathrm{n}}$ and $\mathrm{l}\in {\Omega }_{\mathrm{d}\mathrm{r}\mathrm{p}}$

$$\begin{gathered} {\mathrm{P}}_{\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}= {\mathrm{P}}_{\mathrm{L},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}-{\mathrm{P}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{s},\mathrm{t}}-{\mathrm{P}}_{\mathrm{P}\mathrm{V}}^{\mathrm{p},\mathrm{s},\mathrm{t}}-{\mathrm{P}}_{\mathrm{W}\mathrm{T}}^{\mathsf{\omega },\mathrm{s},\mathrm{t}}+{\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}-{\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{b}\mathrm{a}\mathrm{t}}^{\mathrm{b},\mathrm{s},\mathrm{t}}+{\mathrm{P}}_{\mathrm{c}\mathrm{h},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\mathrm{P}}_{\mathrm{d}\mathrm{i}\mathrm{s},\mathrm{e}\mathrm{v}}^{\mathrm{e},\mathrm{s},\mathrm{t}}+{\mathrm{P}}_{\mathrm{e}\mathrm{l}\mathrm{s}\mathrm{r}}^{\mathrm{h},\mathrm{s},\mathrm{t}}-{\mathrm{P}}_{\mathrm{M}\mathrm{T}}^{\mathrm{h},\mathrm{s},\mathrm{t}}+{\mathrm{P}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{U}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}^{\mathrm{l},\mathrm{s},\mathrm{t}}- {\mathrm{P}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{D}\mathrm{o}\mathrm{w}\mathrm{n}}^{\mathrm{l},\mathrm{s},\mathrm{t}}, \\ \hspace{1em}\hspace{1em}\hspace{1em} \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{v}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{V}\mathrm{p}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}} \end{gathered}$$

(45)

$$\begin{gathered} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathrm{Q}}_{\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}= {\mathrm{Q}}_{\mathrm{L},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}-{\mathrm{Q}}_{\mathrm{D}\mathrm{G}}^{\mathrm{d},\mathrm{s},\mathrm{t}} +{\mathrm{Q}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{U}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}^{\mathrm{l},\mathrm{s},\mathrm{t}}- {\mathrm{Q}}_{\mathrm{L},\mathrm{d}\mathrm{r}\mathrm{p}\mathrm{D}\mathrm{o}\mathrm{w}\mathrm{n}}^{\mathrm{l},\mathrm{s},\mathrm{t}}, \\ \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{v}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{V}\mathrm{p}\mathrm{p}} , \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}} \end{gathered}$$

(46)

$$\left({\mathrm{V}}_{\mathrm{b},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\le {\mathrm{V}}_{\mathrm{b},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}\le \left({\mathrm{V}}_{\mathrm{b},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\right), \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{v}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{V}\mathrm{p}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(47)

$${\left({\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v},\mathrm{s},\mathrm{t}}\right)}^2 +{\left({\mathrm{Q}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v},\mathrm{s},\mathrm{t}}\right)}^2\le {\left({\mathrm{S}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{u}\mathrm{v},\mathrm{M}\mathrm{a}\mathrm{x}}\right)}^2, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{u},\mathrm{v}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{V}\mathrm{p}\mathrm{p}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(48)

3. Proposed Co-Optimization Framework

3.1. Concept Description

The growing adoption of rooftop solar panels, batteries, and other DERs presents both technical and financial opportunities within the power system [8]. This can be enabled by aggregating these resources into VPPs [9]. However, if these VPPs are not handled appropriately by various agents of the power systems, mainly the distribution network service providers, power system infrastructure may run into different technical problems, such as exceeding branch capacity limitations and voltage instability. It is important to note that breaching network security, such as line capacity limits, can lead to significant financial consequences and reduced societal welfare. Previous studies [5,6] have assumed to have direct control over managing the behind-the-meter DERs. However, in countries like Australia [7], DNSPs do not have the authority to access data from behind-the-meter equipment, in order to preserve the privacy of the different layers within the power system.

Co-optimization, which means optimized collaboration among many entities working together to achieve a common goal [27], can be deployed for managing the technical aspects of power layers. In this research, as shown in Figure 1, co-optimization stands for collaboration between DNSPs and VPPs in such a way that makes the operation of the system secure. This will be achieved through the development of a novel co-optimization framework, through which the limits of power trading between such service providers are effectively computed with a guarantee that power trading of VPPs with the markets will not violate the technical aspects of distribution networks. Precisely, under a decentralized framework for preserving the privacy of these mentioned service providers, permissible boundaries are calculated by DNSPs [28].

3.2. Boundary Publishment Development

DNSPs cannot access VPP data in detail because of legislation regarding privacy and independence in operation. Nevertheless, for the security and efficiency of the whole network’s operation, DNSPs should establish and disseminate the operating boundaries of the VPPs, which also requires a framework that will be fast, reliable, and scalable towards changes in dynamic network conditions.

To fill this requirement, a sensitivity analysis framework is developed that will enable DNSPs to calculate and dispatch these boundaries in real-time circumstances. The approach utilizes the relationship between network constraints and the routes of installed VPPs to enable DNSPs to compute boundaries without requiring direct access to VPP data.

As an example of the proposed sensitivity analysis, consider a simplified distribution network integrated with two VPPs, as shown in Figure 2. In this regard, the DNSP is supposed to interact with both VPPs to ensure secure operation of distribution networks without violating their operational limitations. The DNSP must also establish a decentralized and privacy-preserving framework capable of addressing the future needs of energy systems, which increasingly depend on distributed resources, to ensure grid reliability and stability. To achieve this, the following procedures are followed.

Firstly, the DNSP runs the power flow subject to corresponding constraints without the penetration VPPs, which means ${\mathrm{P}}_{\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{t}}=0$ and ${\mathrm{Q}}_{\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{t}}=0$. To accomplish this, the power flow equations of the distribution networks are defined as follows.

$${\mathrm{P}}_{\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n},\mathrm{s},\mathrm{t}}=\sum _{\mathrm{m}\in {\Omega }_{{\mathrm{p}\mathrm{a}\mathrm{r}}_{\mathrm{n}}}}{\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{n},\mathrm{s},\mathrm{t}}-\sum _{\mathrm{k}\in {\Omega }_{{\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{l}}_{\mathrm{n}}}}{\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n}\mathrm{k},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{n}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{n}\mathrm{e}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(49)

$${\mathrm{Q}}_{\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n},\mathrm{s},\mathrm{t}}=\sum _{\mathrm{m}\in {\Omega }_{{\mathrm{p}\mathrm{a}\mathrm{r}}_{\mathrm{n}}}}{\mathrm{Q}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{n},\mathrm{s},\mathrm{t}}-\sum _{\mathrm{k}\in {\Omega }_{{\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{l}}_{\mathrm{n}}}}{\mathrm{Q}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n}\mathrm{k},\mathrm{s},\mathrm{t}},\forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{n}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{n}\mathrm{e}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(50)

$$\begin{gathered} {\mathrm{V}}_{\mathrm{b},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n},\mathrm{s},\mathrm{t}}\approx {\mathrm{V}}_{\mathrm{b},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m},\mathrm{s},\mathrm{t}}-{\displaystyle \frac{\left({\mathrm{R}}_{\mathrm{L},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{n},\mathrm{t}}{\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{n},\mathrm{s},\mathrm{t}}+{\mathrm{X}}_{\mathrm{L}}^{\mathrm{m}\mathrm{n},\mathrm{t}}{\mathrm{Q}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{n},\mathrm{s},\mathrm{t}}\right)}{{\mathrm{V}}_{{\mathrm{s}}_0,\mathrm{n}\mathrm{e}\mathrm{t}}}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{n},\mathrm{m} \\ \hspace{1em}\hspace{1em}\hspace{1em}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{n}\mathrm{e}\mathrm{t}}, {\mathrm{s}}_0 \mathrm{i}\mathrm{s} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{r}\mathrm{o}\mathrm{o}\mathrm{t} \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} \mathrm{o}\mathrm{f} \mathrm{B}\mathrm{u}\mathrm{s} \mathrm{n} \mathrm{a}\mathrm{n}\mathrm{d} \mathrm{m}. \end{gathered}$$

(51)

For bus n which is connected with VPP v, and $\mathrm{v}\in {\Omega }_{\mathrm{V}\mathrm{P}\mathrm{P}}$

$${\mathrm{P}}_{\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n},\mathrm{s},\mathrm{t}}= {\mathrm{P}}_{\mathrm{L},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n},\mathrm{s},\mathrm{t}}+{\mathrm{P}}_{\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{n}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{n}\mathrm{e}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(52)

$${\mathrm{Q}}_{\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n},\mathrm{s},\mathrm{t}}= {\mathrm{Q}}_{\mathrm{L},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n},\mathrm{s},\mathrm{t}}+{\mathrm{Q}}_{\mathrm{V}\mathrm{p}\mathrm{p}}^{\mathrm{v},\mathrm{s},\mathrm{t}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{n}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{n}\mathrm{e}\mathrm{t}} , \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(53)

$$\left({\mathrm{V}}_{\mathrm{b},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\le {\mathrm{V}}_{\mathrm{b},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{n},\mathrm{s},\mathrm{t}}\le \left({\mathrm{V}}_{\mathrm{b},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\right), \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{n}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s}.\mathrm{n}\mathrm{e}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(54)

$${\left({\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{n},\mathrm{s},\mathrm{t}}\right)}^2 +{\left({\mathrm{Q}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{n},\mathrm{s},\mathrm{t}}\right)}^2\le {\left({\mathrm{S}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{m}\mathrm{n},\mathrm{M}\mathrm{a}\mathrm{x}}\right)}^2, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}, \forall \mathrm{n},\mathrm{m}\in {\Omega }_{\mathrm{b}\mathrm{u}\mathrm{s},\mathrm{n}\mathrm{e}\mathrm{t}}, \forall \mathrm{s}\in {\Omega }_{\mathrm{s}\mathrm{e}\mathrm{n}}$$

(55)

After running the above load flow for the distribution networks and calculating the power line conditions, the additional line capacity (ALC) is calculated as follows:

$$\begin{gathered} \mathrm{A}\mathrm{L}\mathrm{C}=\left[\begin{array}{ccc}{∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{\mathrm{1,1}}& \cdots & {∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{1,\mathrm{t}}\\ ⋮& \ddots & ⋮\\ {∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{\mathrm{B}\mathrm{r},1}& \cdots & {∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{\mathrm{B}\mathrm{r},\mathrm{t}}\end{array}\right] \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\mathrm{r}\mathrm{e}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{e}\left(\left[\begin{array}{c}{\mathrm{S}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{1,\mathrm{M}\mathrm{a}\mathrm{x}}\\ ⋮\\ {\mathrm{S}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{B}\mathrm{r},\mathrm{M}\mathrm{a}\mathrm{x}}\end{array}\right],1,\mathrm{t}\right)-\left[\begin{array}{ccc}{\mathrm{S}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{1,1}}& \cdots & {\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{1,\mathrm{t}}\\ ⋮& \ddots & ⋮\\ {\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{B}\mathrm{r},1}& \cdots & {\mathrm{P}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{B}\mathrm{r},\mathrm{t}}\end{array}\right],\forall \mathrm{B}\mathrm{r} \\ \hspace{1em}\hspace{1em}\hspace{1em}\in {\Omega }_{\mathrm{B}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{h}},\forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}} \end{gathered}$$

(56)

In the above expression, $\mathrm{A}\mathrm{L}\mathrm{C}$ means the additional line capacity; ${∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{\mathrm{B}\mathrm{r},\mathrm{t}}$ refers to the available capacity of the branch $\mathrm{B}\mathrm{r}$ at time $\mathrm{t}$ which is calculated through the maximum line capacity minus the line power ${\mathrm{S}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{B}\mathrm{r},\mathrm{t}}$. It is notable that ${\mathrm{S}}_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{n}\mathrm{e}\mathrm{t}}^{\mathrm{B}\mathrm{r},\mathrm{t}}$ is calculated by running the power flow for the distribution networks without considering the VPPs.

The additional line capacity has been calculated so far without the VPP’s penetrations. Now, a sensitivity analysis can be formed to model the behavior of VPPs on the distribution networks. In more detail, a matrix populated with binary numbers (0 and 1) is constructed based on the structure of the distribution network, starting from the substation (SS) and extending to the bus where the VPP is installed.

The sensitivity matrix will be constructed based on a binary system expressing the relationship between branches in the distribution network and VPPs. It will then be identified as follows:

• A branch on the route from the substation to VPP is marked as 1.
• Otherwise, it is 0.

To understand this context, let us go to the modified sample network that has been shown in Figure 3 for illustration. The branches that come along the route of VPP to SS are filled up by 1 and those that do not come under this route are filled by 0. In the example shown in Figure 3, the lines that are filled up with 1 have been highlighted in blue color for easy following.

The highlighted branches, in blue, represent the paths between SS and each VPP.

• In VPP1, branches 1, 2, and 4 give the path.
• In VPP2, flow occurs through branches 1 and 5.

The filled matrix for the mentioned paths is provided as follows.

$$\left[\begin{array}{cc}\stackrel{\mathrm{V}\mathrm{p}\mathrm{p}1}{\overbrace{\begin{array}{c}1\\ 1\\ \begin{array}{c}0\\ \begin{array}{c}1\\ 0\end{array}\end{array}\end{array}}}& \stackrel{\mathrm{V}\mathrm{p}\mathrm{p}2}{\overbrace{\begin{array}{c}1\\ 0\\ \begin{array}{c}0\\ \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}}}\end{array}\right]$$

(57)

• Row i corresponds to branch i.
• Column j corresponds to VPP j.

The above matrix was constructed manually to demonstrate the method. However, to generalize the process and enable it across different networks, a detailed algorithm is provided in Appendix A. This generalized approach systematically identifies the sensitivity matrix, ensuring adaptability to various networks.

By combining Equations (56) and (57), the relationship between the available boundaries of VPPs, the sensitivity analysis matrix, and the additional line capacity can be formulated as a constraint for distribution networks. Specifically, ${\mathrm{B}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}_{\mathrm{V}\mathrm{p}\mathrm{p}1}^{\mathrm{t}}$ and ${\mathrm{B}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}_{\mathrm{V}\mathrm{p}\mathrm{p}2}^{\mathrm{t}}$ are two variables that are optimized to determine the boundaries for VPPs, ensuring that branch capacities are not exceeded.

$$\left[\begin{array}{cc}\stackrel{\mathrm{V}\mathrm{p}\mathrm{p}1}{\overbrace{\begin{array}{c}1\\ 1\\ \begin{array}{c}0\\ \begin{array}{c}1\\ 0\end{array}\end{array}\end{array}}}& \stackrel{\mathrm{V}\mathrm{p}\mathrm{p}2}{\overbrace{\begin{array}{c}1\\ 0\\ \begin{array}{c}0\\ \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}}}\end{array}\right]\left[\begin{array}{c}{\mathrm{B}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}_{\mathrm{V}\mathrm{p}\mathrm{p}1}^{\mathrm{t}}\\ {\mathrm{B}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}_{\mathrm{V}\mathrm{p}\mathrm{p}2}^{\mathrm{t}}\end{array}\right]\le \mathrm{R}\mathrm{e}\mathrm{a}\mathrm{l}\left(\left[\begin{array}{c}{∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{1,\mathrm{t}}\\ ⋮\\ {∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{5,\mathrm{t}}\end{array}\right]\right), \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}$$

(58)

In the above, the explanation was provided for the sample network, but its general form is provided here for different networks.

$${\mathrm{O}\mathrm{F}}_{\mathrm{D}\mathrm{N}}=\mathrm{m}\mathrm{a}\mathrm{x}\left(\sum _{\mathrm{v}\in {\Omega }_{\mathrm{V}\mathrm{p}\mathrm{p}}}\sum _{\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}}{\mathrm{B}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}_{{\mathrm{V}\mathrm{p}\mathrm{p}}_{\mathrm{v}}}^{\mathrm{t}}\right), \forall \mathrm{v}\in {\Omega }_{\mathrm{V}\mathrm{p}\mathrm{p}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}$$

(59)

Subject to

$$\left[\begin{array}{ccc}{\mathrm{R}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{e}}_{\mathrm{V}\mathrm{p}\mathrm{p}1}^1& \cdots & {\mathrm{R}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{e}}_{\mathrm{V}\mathrm{p}\mathrm{p}\mathrm{v}}^1\\ ⋮& \ddots & ⋮\\ {\mathrm{R}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{e}}_{\mathrm{V}\mathrm{p}\mathrm{p}1}^{\mathrm{B}\mathrm{r}}& \cdots & {\mathrm{R}\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{e}}_{\mathrm{V}\mathrm{p}\mathrm{p}\mathrm{v}}^{\mathrm{B}\mathrm{r}}\end{array}\right]\left[\begin{array}{c}{\mathrm{B}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}_{\mathrm{V}\mathrm{p}\mathrm{p}1}^{\mathrm{t}}\\ \begin{array}{c}⋮\\ {\mathrm{B}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}}_{{\mathrm{V}\mathrm{p}\mathrm{p}}_{\mathrm{v}}}^{\mathrm{t}}\end{array}\end{array}\right]\le \mathrm{R}\mathrm{e}\mathrm{a}\mathrm{l}\left(\left[\begin{array}{c}{∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{1,\mathrm{t}}\\ ⋮\\ {∆\mathrm{S}}_{\mathrm{D}\mathrm{N}}^{\mathrm{B}\mathrm{r},\mathrm{t}}\end{array}\right]\right),\forall \mathrm{v}\in {\Omega }_{\mathrm{V}\mathrm{p}\mathrm{p}}, \forall \mathrm{t}\in {\Omega }_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}. \forall \mathrm{B}\mathrm{r}\in {\Omega }_{\mathrm{B}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{h}}$$

(60)

and constraints (49)–(55).

3.3. Procedure of Dynamic Boundary Publishment to the VPPs

Figure 4 illustrates that the process starts by gathering necessary data about the distribution network from the DNSPs, such as the information on network branches and forecasted load demand, among any other relevant operational data. The engine aims at setting the maximum possible boundary for each VPP in such a way that critical network constraints are satisfied.

Once the optimization engine is run, it provides the available capacities for each VPP, indicating the boundary each VPP can be allocated without violating network constraints. This value is then sent back to the VPPs through a secure data exchange hub. In turn, the VPPs take the assigned boundary limits to schedule their resources in such a way to achieve maximum profit within the given boundary.

After scheduling, the VPPs send their bids to the DNSPs. The DNSPs will re-run the load flow with the inclusion of bids coming from the VPPs to check the technical aspects of the network are not breached. This makes sure that, in summary, the overall network constraints remain within the standard limits and keep the grid stable.

The pseudocode of the proposed method is provided in

Table 1. The pseudocode of the proposed boundary publishment approach.

The Pseudocode of Dynamic Boundary Publishment.
Step 1:	Get data related to the distribution networks, such as loads, configuration, line limitations, etc.
Step 2:	Run load flow for the distribution network and calculate the additional line capacity by Equation (56).
Step 3:	Construct the sensitivity matrix based on the lines that connect the VPPs to the distribution network based on the algorithm provided in the Appendix A.
Step 4:	Maximize the boundaries of VPPs by Equation (59) subject to the constraints described in Equation (60).
Step 5:	DNSP Publishes the boundaries to the VPPs via data exchange hub.
Step 6:	The service providers of VPPs receive the boundaries to schedule their DERs with the aim of maximizing profit.
Step 7:	DNSPs receive bids from VPPs and perform power flow analyses in the presence of these bids. This ensures that the system’s technical constraints are not breached while minimizing voltage deviations or operational costs.
Step 8:	End.

below to outline the steps involved.

4. Results and Discussions

This section presents the simulation results demonstrating the performance of the proposed voltage and congestion management approach when the DNSP publishes the operational boundaries to the VPPs. The methodology is applied to the IEEE 33-bus radial distribution network [29], which is connected to two VPPs, as shown in Figure 5. The corresponding demands are illustrated in Figure 6.

Two VPPs are formed based on a 12-bus IEEE network [30] and each consists of various resources, including renewable energy sources, hydrogen fuel cell demands, market prices, and load factors, all of which are presented in Figure 7. Each VPP includes dispatchable generators (DGs) and microturbines (MTs) with capacities of 200 kW and 100 kW, respectively. Additional details regarding the storage facilities and relevant price coefficients can be found in

Table 2. Battery storage system characteristics.

Capacity	${\mathit{P}}_{\mathit{c}\mathit{h},\mathit{b}\mathit{a}\mathit{t}}^{\mathit{b},\mathit{m}\mathit{a}\mathit{x}}$	${\mathit{P}}_{\mathit{d}\mathit{i}\mathit{s},\mathit{b}\mathit{a}\mathit{t}}^{\mathit{b},\mathit{m}\mathit{a}\mathit{x}}$	${\mathit{\eta }}_{\mathit{c}\mathit{h},\mathit{b}\mathit{a}\mathit{t}}$	${\mathit{\eta }}_{\mathit{d}\mathit{i}\mathit{s},\mathit{b}\mathit{a}\mathit{t}}$
4 MWh	0.4 MW	0.4 MW	0.95	0.96

,

Table 3. EVA characteristics.

Capacity	${\mathit{P}}_{\mathit{c}\mathit{h},\mathit{e}\mathit{v}}^{\mathit{e},\mathit{m}\mathit{a}\mathit{x}}$	${\mathit{P}}_{\mathit{d}\mathit{i}\mathit{s},\mathit{e}\mathit{v}}^{\mathit{e},\mathit{m}\mathit{a}\mathit{x}}$	${\mathit{P}}_{\mathit{t}\mathit{r}\mathit{i}\mathit{p},\mathit{e}\mathit{v}}^{\mathit{e},\mathit{t}}$	${\mathit{\eta }}_{\mathit{c}\mathit{h},\mathit{e}\mathit{v}}$	${\mathit{\eta }}_{\mathit{d}\mathit{i}\mathit{s},\mathit{e}\mathit{v}}$	Trip Times
2 MWh	0.2 MW	0.2 MW	0.02 MW	0.95	0.96	7, 16, 24

,

Table 4. Hydrogen system characteristics.

Hydrogen Storage Capacity	${\mathit{H}}_{\mathit{e}\mathit{l}\mathit{s}\mathit{r}}^{\mathit{h},\mathit{m}\mathit{a}\mathit{x}}$	${\mathit{H}}_{\mathit{d}\mathit{i}\mathit{s},\mathit{h}\mathit{g}\mathit{n}}^{\mathit{h},\mathit{m}\mathit{a}\mathit{x}}$	${\mathit{\eta }}_{\mathit{e}\mathit{l}\mathit{s}\mathit{r}}$	${\mathsf{\Psi }}_{\mathit{e}\mathit{i}\mathit{h}}$
100 kg	10 kg	10 kg	0.8	$39\mathrm{k}\mathrm{W}\mathrm{h}/\mathrm{k}\mathrm{g}\mathrm{H}2$

and

Table 5. The cost coefficients for different DERs.

${\mathit{C}}_{\mathit{D}\mathit{G}}^{\mathit{d}}$	${\mathit{C}}_{\mathit{O}\&\mathit{M},\mathit{M}\mathit{T}}^{\mathit{h}}$	${\mathit{C}}_{\mathit{c}\mathit{y}\mathit{c},\mathit{b}\mathit{a}\mathit{t}}^{\mathit{b}}$	${\mathbf{C}}_{\mathit{r}\mathit{e}\mathit{m}}^{\mathit{e}}$	${\mathit{C}}_{\mathit{f}\mathit{c}\mathit{v}}^{\mathit{h}}$	${\mathit{C}}_{\mathit{O}\&\mathit{M},\mathit{e}\mathit{l}\mathit{s}\mathit{r}}^{\mathit{h}}$	${\mathit{C}}_{\mathit{O}\&\mathit{M},\mathit{P}\mathit{V}}^{\mathit{p}}$	${\mathit{C}}_{\mathit{O}\&\mathit{M},\mathit{W}\mathit{T}}^{\mathit{\omega }}$	${\mathbf{C}}_{\mathit{c}\mathit{m}\mathit{r}}^{\mathit{c}}$
0.06 ($/kWh)	0.015 ($/kWh)	0.30 ($/switching)	0.005 ($/kWh)	9 ($/kgH2)	1.50 ($/kgH2)	0.005 ($/kWh)	0.01 ($/kWh)	0.12 ($/kWh)

. Two cases are considered to evaluate the effectiveness of the proposed approach:

• Case 1: No boundaries are set for the VPPs, allowing unrestricted power trading between the VPPs and the distribution network.
• Case 2: Boundaries are issued by the DNSP, limiting power trading between the VPPs and the distribution network.

4.1. Results from DNSP

In this section, the voltage levels and power flow through the branches of the distribution network are analyzed for both Case 1 and Case 2. Specifically, in Case 1, the power trading of the VPPs is unrestricted within the distribution network. In contrast, Case 2 operates under the boundaries published by the DNSP, based on the proposed method, as shown in Figure 8 for both VPP1 and VPP2.

From a voltage management perspective, Figure 9 presents the voltage levels of the distribution network as box plots for both cases. In Case 1, voltage breaches the standard limit of 0.95 p.u. from hours 10 to 18 due to unregulated power trading between the DNSP and VPPs. These voltage fluctuations are severe, risking damage to network infrastructure and reducing the reliability of the network. However, in Case 2, where the published limits set by the DNSP are followed (see Figure 8), voltages remain within the permissible range of 0.95–1.05 p.u. It is observed that there is greater variation in the voltage between buses 18 and 19, as bus 18 is farther from the substation bus while bus 19 is closer, leading to more noticeable fluctuations.

In terms of congestion management, Figure 10 shows the power flows in network branches as box plots for both cases. In Case 1, some branches exceed their capacity limits due to unconstrained power trading with the VPPs. For example, from branch 6 to 17, the maximum capacity limits are breached, posing significant technical risks such as overheating and potential equipment damage. This highlights the critical role of congestion management in preventing such issues. Conversely, in Case 2, where the VPPs adhere to the boundaries set by the DNSP, no congestion occurs and the power flow stays within acceptable limits.

The comparative analysis underscores the importance of a co-optimization approach between DNSPs and VPPs. Without coordination, the system is vulnerable to voltage fluctuations and congestion within distribution networks. This approach not only ensures reliable power delivery but also aligns with regulatory frameworks designed to maintain system integrity in increasingly complex and decentralized energy markets.

In addition to technical assessments, this study includes an economic evaluation to highlight how the profit of VPPs is affected when adhering to the boundaries set by the DNSP. To illustrate this, Figure 11 presents the profits of VPPs under different cases. As shown, the profit of VPP₁ decreases by approximately 20% when constrained by the DNSP’s boundaries. In contrast, VPP₂ experiences a smaller reduction of about 7%. This disparity is primarily due to the stricter boundaries imposed on VPP₁, as illustrated in Figure 8. In other words, as shown in Figure 5, VPP₁ is connected to bus 18, which is further away from the substation. This location results in voltage levels closer to the lower limit (0.95 p.u.), necessitating tighter boundaries to maintain compliance with voltage standards. On the other hand, VPP₂, connected to bus 21, is located closer to the substation. Its higher voltage levels allow for more permissive boundaries. Overall, the findings indicate that the profits of VPPs are not significantly impacted when adhering to the DNSP-provided boundaries, ensuring the safe operation of distribution networks.

4.2. Results from VPPs

From the VPP perspective, the voltage profiles of internal networks for both VPPs are shown in Figure 12 and Figure 13, respectively. From these figures, it is observed that the voltage levels of the VPPs are within the range of acceptability. However, there are noticeable differences in voltage patterns between the two cases. As can be seen, in Case 1, for intervals 1 to 16, the voltage of the VPPs is less than in Case 2. At hour 17 up to 24, however, the voltage level of the VPPs is higher in Case 1 than that in Case 2. This is because, in Case 1, the absence of boundaries for the VPPs enables them to buy power either from the market or excess generation of renewables for charging their storage facilities at hours 1–16 and sell between hours 17 and 24, which leads to increased voltage. However, for Case 2, with dynamic boundaries limiting power trading, the VPPs must rely more on their local DERs. Because of the storage facilities within the VPPs, they have flexibility in managing their resources effectively and maintaining the technical aspects of the VPPs.

When analyzing the scheduling of the VPPs’ resources, Figure 14 and Figure 15 depict various patterns. In Case 1, both VPPs purchase significant power from the markets during hours 2 and 16, while selling during intervals 17 to 24 when market prices are higher to maximize profit. In contrast, in Case 2, these power trading quantities decrease as they must comply with the limits imposed by the DNSP, leading to an increased reliance on the local resources to meet the demand. For example, in Case 2, the DG of VPP1, represented by the grey bar, plays an important role in load supply. The same happens for VPP2, which is also constrained to be within the published boundaries by the DNSP for Case 2. Other resources are used to maintain the power balance for the different periods. The main difference here is that the boundaries for VPP2 are higher; hence, it can trade more power with the market and has a higher profit potential.

If we examine the hydrogen system scheduling in Figure 16 and Figure 17, it becomes evident that the electrolyzer possesses different generation patterns in the two cases. In Case 1, the electrolyzer generates hydrogen continuously. Conversely, in Case 2, significant hydrogen generation occurs during intervals with PV generation, while generation during other intervals is reduced because of having the boundaries on power trading with the distribution network. As a result, the VPP schedules its hydrogen system based on local resources to maximize profit, leveraging power-to-hydrogen to capture surplus PV generation and supplying the hydrogen fuel cell vehicles.

4.3. Discussion and Summary

A brief explanation of the assumptions of the mathematical model are presented here and they are evaluated for their application in a real-world context. Since the complexity of modern electricity markets is growing, fast methods are needed to for real-time market applications. One of the most important factors in assessing the proposed methodology is computational efficiency. The time consumption of the solver was measured on an Intel core i7 laptop with 16 GB of memory, using GUROBI version 11 as the solver. More importantly, it only takes 0.5 s to publish the boundaries to the VPPs, illustrating well the merit of relying on convex programming and serving present electricity market requirements. In particular, compared to population-based algorithms (such as genetic algorithms), convex programming brings some crucial advantages. First, it is much faster. The genetic algorithms are iterative and take many iterations for convergence, which is unacceptable for a real-time market. Genetic algorithms have no guarantees of reaching the global solution and hence might end up finding a local solution. In a nutshell, the proposed method outperforms population-based algorithms in both the time efficiency and accuracy of the solution; hence, convex optimization is more suitable for real-time electricity market operation.

According to the obtained results, the key findings of this study are outlined in

Table 6. The key findings of the investigation.

Comprehensive VPP Model	Integrates various DERs, including hydrogen and electricity, for increased flexibility.
Decentralized operation	Ensures privacy and confidentiality of service providers’ data, in line with emerging market privacy requirements.
Network safety assurance	It mitigates the risk of branch congestion by publishing boundaries to VPPs.
Voltage compliance	Keeps voltage above 0.95 p.u., guaranteeing compliance with operational standards after the method is implemented.
Real-time boundary publication	Support for real-time boundary publishing allows integration into real-time markets because of the fast response in just 0.5 s.

, which highlights the significant outcomes of implementing the proposed co-optimization method.

5. Conclusions

This paper introduces a co-optimization technique that enables DNSPs to dynamically set and publish the operational boundaries for VPPs while ensuring effective management of their interactions, without compromising the integrity of the distribution network. Operating in a decentralized manner, DNSPs determine these boundaries without disclosing any private VPP data, thereby maintaining the confidentiality of all service providers. The proposed approach offers several key advantages:

• The VPP model accommodates a wide range of distributed energy resources, including hydrogen and electricity, which enhances the flexibility and adaptability of the VPP.
• The decentralized nature of the method ensures that service provider privacy is upheld, aligning with the emerging regulatory requirements of modern energy markets.
• The co-optimization framework guarantees safe network operation by mitigating risks such as voltage instability and branch congestion. For example, in the unbounded case study, the network voltage fell below the minimum acceptable threshold of 0.95 p.u. However, after implementing the proposed approach, voltage levels were maintained above 0.95 p.u., ensuring compliance with operational standards.
• The convex nature of the proposed framework improves computational efficiency, making it compatible with advanced solvers. The non-iterative structure enables real-time boundary publication, supporting seamless integration into online market operations. The solver requires only 0.5 s to publish boundaries using a standard personal laptop.

While this research provides valuable contributions to voltage and congestion management, it does not address the development of a market-based framework. Specifically, it does not consider VPPs as price makers in market-clearing processes, which presents an opportunity for future research. Additionally, the integration of machine learning techniques to further enhance the boundary publication process, particularly in dynamic environments such as rapid load changes, is another promising avenue for future work