a Abstract-This paper presents a novel measurement-based frequency regulation scheme that utilizes the contributions of inverter-based resources (IBRs). IBRs are assumed to be aggregated within dynamic virtual power plant (DVPP), conceptualized as a new entity in modern power systems that can be summoned by the system operator to provide dynamic ancillary services. While the participation of DVPPs in frequency regulation is outlined in ancillary service markets, the main objective is to distribute regulation signals among the resources within DVPPs. To achieve this, a power imbalance propagation paradigm is introduced, aiming to estimate the extent of power imbalance across the grid. Subsequently, the contribution factor is established for the purpose of evaluating the involvement of individual clusters of IBRs within the context of both primary and secondary frequency control commitments. The proposed mathematics-based approach is scalable and only needs to have access to few local measurments and exchanged power information between the adjacent neighboring areas to realize a real-time application. The effectiveness of the proposed methodologies is examined on two power grids.

I. INTRODUCTION

requency regulation mechanisms play an important role in ensuring reliable and secure operation of electric power systems. As a normal part of conducting business, vertically integrated utilities i.e., large-scale synchronous generators (SGs), had the responsibility to maintain a near real-time balance between generation and load. Following deregulation in the power industry, the responsibility for frequency regulation has been redefined as an ancillary service, which is performed by parties that generate, control, and transmit electricity in support of the basic services of energy supply and power delivery [START_REF] Kirby | Frequency regulation basics and trends[END_REF][START_REF] Björk | Dynamic virtual power plant design for fast frequency reserves: Coordinating hydro and wind[END_REF][START_REF] Bevrani | Grid Connected Converters: Modeling, Stability and Control[END_REF]. Transition towards renewable energy sources (RESs), nevertheless, imposes that regulatory services in future power systems should be shouldered by a combination of SGs and inverter-based resources (IBRs).

Several investigations have been reported in recent years, pertaining to the application of IBRs in frequency regulation. The research conducted so far in the field can be broadly categorized into two groups: 1) model-oriented approaches and 2) model-free-based approaches. Model-oriented methods primarily focus on studying the participation of a limited Enhanced Frequency Regulation Scheme: An Online Paradigm for Dynamic Virtual Power Plant Integration Hêmin Golpîra, Senior Member, IEEE, Bogdan Marinescu, Member, IEEE F number of individual IBRs [START_REF] Ademola-Idowu | Frequency stability using MPC-based inverter power control in low-inertia power systems[END_REF][START_REF] Stanojev | MPC-based fast frequency control of voltage source converters in low-inertia power systems[END_REF][START_REF] Kumar | Decentralized reinforce-control of renewable dynamic virtual power plant enabling it to offer ancillary services: An attempt towards net-zero targets[END_REF][START_REF] Häberle | Control design of dynamic virtual power plants: An adaptive divide-and-conquer approach[END_REF], or aggregated groups of small-scale RESs envisioned within the concept of virtual power plants (VPPs) [START_REF] Dall'anese | Optimal regulation of virtual power plants[END_REF][START_REF] Anderson | Frequency regulation with heterogeneous energy resources: A realization using distributed control[END_REF][START_REF] Espinosa | Referencetracking control policies for packetized coordination of heterogeneous DER populations[END_REF][START_REF] Gu | An Auto-Tuned Robust Dispatch Strategy for Virtual Power Plants to Provide Multi-Stage Real-time Balancing Service[END_REF], within the context of frequency regulation. The increasing integration of RESs, however, emphasizes the significance of regulation methods that efficiently handle operational uncertainties, a capability lacking in model-oriented approaches. More precisely, modeloriented controllers present significant real-time computational obstacles. Furthermore, frequency regulation policies necessitate the development of appropriate disaggregation-based approaches for VPPs to accurately ascertain the contribution of each source in adhering to the instantaneous regulation request, an aspect that is overlooked in recent research. Other shortcomings of the recent modeloriented studies encompass overlooking the dynamic of the host grid, disregarding the interplay between VPPs and grid frequency regulation services, and omitting consideration of RESs at the transmission level.

Model-free category mainly focuses on using learning-based techniques to determine participation of IBRs in the frequency regulation service. The utilization of IBRs for frequency regulation is addressed through reinforcement learning (RL)based approaches in [START_REF] Yan | Data-driven load frequency control for stochastic power systems: A deep reinforcement learning method with continuous action search[END_REF][START_REF] Sadeghi | Optimal bidding strategy of a virtual power plant in day-ahead energy and frequency regulation markets: A deep learning-based approach[END_REF]. A two-stage deep RL-based approach is discussed in [START_REF] Yi | An Improved Two-Stage Deep Reinforcement Learning Approach for Regulation Service Disaggregation in a Virtual Power Plant[END_REF] for predicting the output power of sources and issuing disaggregated request balancing signals. The primary difficulty when employing RL lies in ensuring that the learned controllers are stabilizing, resulting in frequency within the system eventually reaching a stable equilibrium following disturbances. For this purpose, existing approaches often adopt soft penalties, imposing substantial costs when states deviate from predefined ranges. Nevertheless, it is crucial to treat stability as a hard constraint rather than relying on penalties, as training can only be conducted with a limited number of samples, while the controller must ensure stability across a broader set of points in the state space. Furthermore, the typical practice involves using generated trajectories to train the neural network controllers, yet the evolution of state variables over lengthy time horizons makes direct back propagation inefficient. Moreover, challenges associated with algorithmic stability, which concerns how slight variations in input data can impact the output of a machine learning-based approach, restrict the feasibility of employing RL-based methods in real-time applications.

Lyapunov function is incorporated with RL in [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] to guarantee frequency stability. However, simply thresholding the output once the constraints are reached tend to lead to very poor performances [START_REF] Polycarpou | On-line approximation control of uncertain nonlinear systems: issues with control input saturation[END_REF]. More precisely, a primary challenge associated with this method arises when the control signal generated by the control law cannot be fully executed due to input saturation, imposing constraints on the magnitude of the control input. Furthermore, this method tries to maintain the nodal frequency of buses within a desired safe region and to recover from undesired initial conditions while could not guarantee secure operation in transient period. The need for massive datasets to train on, along with the complexity of implementation, are other important drawbacks. For more on VPP and IBRs in the context of ancillary services see [START_REF] Ochoa | Control Systems for Low-Inertia Power Grids: A Survey on Virtual Power Plants[END_REF][START_REF] Yan | Frequency Control and Optimal Operation of Low-Inertia Power Systems with HVDC and Renewable Energy: A Review[END_REF]. This paper addresses the aforementioned drawbacks by proposing a mathematical-based framework that offers frequency regulation for the power grid through appropriately dispatching the resources within a dynamic VPP (DVPP). The DVPP, which was initially conceptualized in the POSYTYF project [START_REF] Marinescu | POSYTYF concept and objectives[END_REF][START_REF] Marinescu | Dynamic virtual power plant: A new concept for grid integration of renewable energy sources[END_REF], is assumed to comprise a diverse range of resources, including conventional SGs and IBRs across different power levels and time scales. The proposed modelfree based method relies on the fundamental mathematical equations governing power systems to generate reference balancing signals for the resources within a DVPP through primary frequency control (PFC) and secondary frequency control (SFC) loops.

To achieve the aforementioned objectives, a measurementbased approach is proposed to determine the power imbalance size at the connecting point (CP) of the resources within the DVPP to the main grid. The participation of resources within the DVPP in meeting the established frequency response reserve (FRR) is subsequently calculated based on the RoCoF and steady-state frequency deviation. The next step involves adjusting the reference power of each resource based on the steady-state frequency deviation through SFC. Finally, an algorithm is proposed for real-time updating of the disaggregated reference balancing signals in response to changing operating points. Therefore, the main contributions of this paper are as follows:

• Introducing the power imbalance propagation paradigm, a measurement-based approach to mathematically model the propagation of power imbalance in the grid. • Proposing a mathematical-based approach to allocate PFC and SFC between the IBRs within a DVPP and adjust the setpoints based on real-time measurements. It can finely dispatch the resources, ensuring secure operation by maintaining the frequency dynamics, including the rate of change of frequency (RoCoF) and steady-state deviation, within a safe range. The approach is scalable and requires access to a limited number of local measurements and exchanged power data within contiguous neighboring areas, in order to facilitate the execution of a real-time application.

II. PROPOSED FREQUENCY REGULATION METHODOLOGY

A. Problem Definition Consider a power system which is subjected to a physical disturbance, causing power imbalance in the system. The system operator (SO) is going to regulate the frequency by utilizing the control loops of both SGs and DVPPs, including PFC and SFC loops, as schematically represented in Fig. 1. The model in Fig. 1 is based on the notation from [START_REF] Lasseter | Grid-forming inverters: A critical asset for the power grid[END_REF], where it is developed as a function of IBRs penetration level denoted as "r". To maintain a consistent 1 pu system load, the power rating of the aggregated model of SGs in the figure decreases by (1-r). It is noteworthy that the simplified model in Fig. 1 is solely for analysis purposes, and all simulations use a 6 th -order generator model with turbine-governing systems for SGs, whereas dynamic models are applied to IBRs.

With reference to the notation in Fig. 1, one could describe the distinct participation of DVPP in both PFC and SFC as

∆𝑃 𝐷𝑉𝑃𝑃 = ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 + ∆𝑃 𝐷𝑉𝑃𝑃 𝑆𝐹𝐶 (1𝑎)
and It is important to highlight that the distinct functionalities of PFC and SFC imply that during PFC, only the terms associated with it on the right-hand side of (1) have non-zero values, while both terms remain non-zero during SFC. PFC: The PFC is responsible for contributing to arresting frequency deviations through the provision of frequency response. Conventionally, the frequency response is defined as the sum of the inertial response, natural load response, and governor response of SGs to frequency deviation, as equivalently illustrated by HG, DG and RG in Fig. 1. The commonly employed strategy for utilizing IBRs to provide PFC involves configuring them to mimic SGs, employing emulated inertial response (HDVPP), emulated load response (DDVPP), and f/p droop control, as illustrated in Fig. 1. However, conventional PFC mechanisms in SGs operate as autonomous loops that don't rely on external inputs to be triggered, whereas inverters are solid-state electronic devices capable of tracking power setpoints. Hence, IBRs within DVPP have the capability to rapidly respond to control signal, denoted as ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 in (1b), promising the implementation of real-time approaches for PFC. This suggests that it is feasible to implement alternative control laws to ensure the stability of an IBR-integrated power system, surpassing a linear droop response, as proposed in other studies such as [START_REF] Stanojev | MPC-based fast frequency control of voltage source converters in low-inertia power systems[END_REF][START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF][START_REF] Johnson | Synchronization of parallel single-phase inverters with virtual oscillator control[END_REF]. This possibility is investigated in this paper, which introduces a new control law, i.e., u(ꞷpll), as a replacement for the current linear f/p control law.

∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 = ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 + ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 (1𝑏)
The SO conducts contingency analysis within specified time intervals, ranging from a few minutes to several hours, with the aim of communicating the permissible frequency and the minimum necessary quantity of FRR, in terms of MW, to the ancillary service market [START_REF] Kontis | Provision of primary frequency response as ancillary service from active distribution networks to the transmission system[END_REF]. From the SO's point of view, DVPP is responsible for providing a fraction of FRR, denoted as ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 , specified based on the ancillary service market settlement. From DVPP owner's point of view, however, the resources within their territory should be dispatched appropriately to fulfil ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 requirement. To achieve this, ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 should be represented based on the emulated frequency response characteristics. Furthermore, DVPP's contribution to inertial response plays a pivotal role in the allocation of resources within DVPP. The contribution of each cluster of IBRs within the DVPP in achieving aggregated frequency response characteristics or inertia constant determines setpoint signals for each resource.

SFC: Once the PFC has arrested and stabilized frequency decline following a disturbance, SFC comes into play to restore the frequency to its nominal value using automatic generation control (AGC) system. The AGC system utilizes system frequency and net actual interchange deviations to determine the area control error (ACE) signal, which indicates whether a system is in balance or requires adjustments to generation resources. The AGC system sends real-time reference signal ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 that prompts the DVPP to perform SFC and counteract the ACE. Again, from DVPP owner's point of view, the reference signal needs to be disaggregated among the resources within the DVPP.

B. Definitions and Hypothesis

This paper examines various hypotheses and definitions to formulate the proposed approach for frequency regulation. These encompass: 1-While a VPP exclusively handles the power injection and integration of individual RES generators into the market, the DVPP concept facilitates the integration of RES generators into the electrical grid by leveraging their collective flexibility. This involves swiftly adjusting output for frequency control, internally managing fluctuations, and marketing their combined generation output in the wholesale market. Compared to the conventional VPP, the DVPP concept will jointly address: ➢ Static power dispatch suggests that the DVPP has the capability to manage the dispatch of all IBRs functioning in grid-forming mode, ensuring optimal operation; ➢ Dynamic power control establishes dynamic regulation for each IBR, as well as for the DVPP, ensuring participation in grid ancillary services. 2-The DVPP is essentially defined as a group of generating units without specific criteria for its formation. 3-The DVPP includes non-dispatchable RESs operating in grid-following mode, and dispatchable ones functioning in grid-forming mode, realized using the proposed control laws. Nonetheless, the DVPP appears to the grid as a cohesive, dispatchable, and fully manageable energy source that can be fine-tuned to support ancillary services, as shown in Fig. 1. The RESs employed in this paper consist of photovoltaics (PVs) and wind generation. 4-The DVPP transmits PFC regulation instructions to IBRs every second in response to frequency fluctuations. To fully implement the SFC loop within the DVPP, it transmits SFC regulation instructions to IBRs at two distinct time intervals. The first action operates on a timescale of 4 to 5 minutes, ensuring optimal production levels of IBRs within the DVPP. The second component involves an internal power redispatch mechanism, operating at an intermediate time-scale, approximately 4 seconds, and is the focus of this paper [START_REF] Adabi | Direct Participation of Dynamic Virtual Power Plants in Secondary Frequency Control[END_REF]. These update periods account for the delay linked to transmitting instructions through cellular wireless communications. 5-The DVPP does not take into account dispatchable loads.

C. Proposed method structure

The proposed frequency regulation framework consists of three main modules: (1) estimating the power imbalance size in the CP of resources to the main grid, called power imbalance propagation paradigm; (2) generating PFC and AGC reference balancing signals and disaggregating the signals between the resources within DVPP; and (3) real-time re-adjusting of the setpoint signals. Fig. 2 illustrates a schematic representation of the proposed frequency regulation, with emphasis on the key modules of interest highlighted in orange, providing a breakdown of the input and output signals for each module.

D. Power Imbalance Propagation Paradigm

Let to consider a power grid consisting n cluster of IBRs connected to the main grid at different voltage levels, as shown in Fig. 3, creating a DVPP structure. It is assumed that clusters 1, 2, …, m are at the distribution level, connected to the main grid via transformers T1,…, Tm, and clusters m+1, …, n are directly connected to the transmission level. The occurred power imbalance in the system, denoted as ∆𝑃 𝐿 , may be manifested at the CP of cluster i, connected to the grid either through transformer Ti or directly in Fig. 3, to the grid with a magnitude of ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 . The preliminary aim is to derive a relationship that can be used to calculate ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 based on a few locally measurable quantities.

The rationale behind developing a framework for power imbalance propagation is rooted in the concept of local frequency estimation in modern power systems. This suggests that frequency is no longer considered a universal property of the grid and can vary from one bus to another [START_REF] Ghosh | An Analytical Approach for Frequency Estimation of Modern Power Grid[END_REF][START_REF] Milano | Frequency divider[END_REF][START_REF] Golpîra | A center-of-gravity-based approach to estimate slow power and frequency variations[END_REF]. Accordingly, one can express the local emulated motion of the under-study cluster using the swing equation of the form [START_REF] Golpîra | Renewable Integrated Power System Stability and Control[END_REF]

2𝐻 𝐷𝑉𝑃𝑃 𝑖 𝑑 2 𝛿 𝑖 𝑑𝑡 2 = 𝑇 𝑚 𝑖 (𝑡) -𝑇 𝑒 𝑖 (𝑡) ⏟ ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 (2)
where 𝛿 𝑖 , 𝐻 𝐷𝑉𝑃𝑃 𝑖 , 𝑇 𝑚 𝑖 and 𝑇 𝑒 𝑖 express the emulated quantities of rotor angle, inertia constant, input mechanical torque, and output electrical torque, respectively, for cluster i. Furthermore, ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 refers to the size of power imbalance at CP of cluster i to the main grid. By reorganizing the terms in (2) based on the defined expression for RoCoF as 𝑅𝑜𝐶𝑜𝐹 = 𝑑(𝛥𝜔) 𝑑𝑡 , one obtains

2𝐻 𝐷𝑉𝑃𝑃 𝑖 × 𝑅𝑜𝐶𝑜𝐹 𝑖 = ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 (3)
Since ( 3) is derived based on the concept of a virtual SG for cluster i, the concepts of frequency and RoCoF are defined as the primary and secondary derivative of the measured voltage phase angle, respectively. Hence, by calculating the contribution of cluster i to the equivalent inertia constant of the grid, denoted as 𝐻 𝐷𝑉𝑃𝑃 𝑖 , one could estimate ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 .

To calculate 𝐻 𝐷𝑉𝑃𝑃 𝑖 , let us define the injected power of cluster i, operated in the grid-forming mode, to the main grid as 𝑃 𝑖 𝐺𝑟𝑖𝑑 . By applying 𝑓 𝑖 -𝑃 𝑖 𝐺𝑟𝑖𝑑 data to the Prony analysis, one obtains

𝑓 𝑖 𝑃 𝑖 𝐺𝑟𝑖𝑑 = ∑ 𝑎 𝑗 𝑠 + 𝑐 𝑗 𝐾 𝑗=1 = ∑ 1 2𝐻 𝑗 𝑖 𝑠 + 𝐷 𝑗 𝑖 𝐾 𝑗=1 (4)
where 𝑓 𝑖 , 𝑃 𝑖 𝐺𝑟𝑖𝑑 , 𝐻 𝑗 𝑖 and 𝐷 𝑗 𝑖 are frequency at the CP to the grid, injected power to the main grid, components of emulated inertia constant for cluster i as seen from the grid, and components of emulated damping factor for cluster i as seen from the grid, respectively. Eq. ( 4) reveals that K equivalent SGs with characteristics represented by 𝐻 𝑗 𝑖 can emulate the aggregated dynamic behavior of IBRs in cluster i, as viewed from the perspective of the main grid. Since these equivalent SGs are assumed to be located on the same bus as the under-studied cluster of IBRs, they exhibit complete coherence. Consequently, the emulated inertia constant of the cluster can be represented by the equivalent inertia for a group of K coherent SGs, namely

𝐻 𝐷𝑉𝑃𝑃 𝑖 = ∑ 𝐻 𝑗 𝑖 𝐾 𝑗=1 (5)
It is essential to emphasize that the analytical equation in ( 5) is based on locally measured quantities instead of relying on mechanical momentum, making it readily applicable to either SGs or clusters of IBRs.

Equations ( 3) through ( 5) could be interpreted as a fault propagation paradigm which could be utilized to estimate the size of power imbalance throughout the grid. Although [START_REF] Bevrani | Grid Connected Converters: Modeling, Stability and Control[END_REF] appears to decouple the dynamics of cluster i from the rest of the grid, the impact of the grid and disturbances is actually manifested in (4), [START_REF] Stanojev | MPC-based fast frequency control of voltage source converters in low-inertia power systems[END_REF], and local RoCoF measurements.

The subsequent step involves allocating the balancing power associated with PFC and SFC to the individual resources within the DVPP.

E. Frequency Regulation Framework

1) PFC Allocating PFC commitment between the resources within the DVPP should be approached from two distinct perspectives: 1) a dynamic point of view which determines the supplementary contribution of DVPP in inertial response of the system and 2) a steady-state point of view which determines supplementary contribution of DVPP in frequency response characteristic of the system.

Dynamic point of view: The RoCoF at the center of inertia (COI) of the SG-based power grid, denoted as 𝑅𝑜𝐶𝑜𝐹 𝐶𝑂𝐼 , after the occurrence of the power imbalance ∆𝑃 𝐿 can be computed using (3) as

𝑅𝑜𝐶𝑜𝐹 𝐶𝑂𝐼 = ∆𝑃 𝐿 2𝐻 𝐺 (6)
where 𝐻 𝐺 is the equivalent grid inertia constant.

Given the involvement of DVPP in the PFC, one could express the 𝑅𝑜𝐶𝑜𝐹 𝐶𝑂𝐼 based on the equivalent inertia 𝐻 𝐺 + 𝐻 𝐷𝑉𝑃𝑃 , as

𝑅𝑜𝐶𝑜𝐹 𝐶𝑂𝐼 = ∆𝑃 𝐿 2(𝐻 𝐺 + 𝐻 𝐷𝑉𝑃𝑃 ) (7)
Secure operation of the system requires that the RoCoF following the occurrence of the power imbalance should remain within a pre-specified range, as defined in the grid codes, namely where 𝑅𝑜𝐶𝑜𝐹 𝑑𝑒𝑠 is the desired RoCoF. Re-arranging the term in (8) yields to

𝐻 𝐷𝑉𝑃𝑃 ≥ ∆𝑃 𝐿 2𝑅𝑜𝐶𝑜𝐹 𝑑𝑒𝑠 -𝐻 𝐺 (9)
Indeed, the right-hand side of ( 9) specifies the minimum contribution of DVPP in the system to ensure secure operation, adhering to the RoCoF requirements. Here, it should be noted that IBRs operating in a grid-forming mode may enhance the system's overall inertia using the virtual SG concept, making it feasible to adjust the inertia level for a specific system in real-time, as addressed in research works such as [START_REF] Golpîra | Emulation of virtual inertia to accommodate higher penetration levels of distributed generation in power grids[END_REF][START_REF] Kryonidis | Use of ultracapacitor for provision of inertial response in virtual synchronous generator: Design and experimental validation[END_REF]. Consequently, equation ( 9) gains practical significance.

To calculate the contribution of each cluster of resources within the DVPP to realize 𝐻 𝐷𝑉𝑃𝑃 [START_REF] Anderson | Frequency regulation with heterogeneous energy resources: A realization using distributed control[END_REF], consider the i th cluster in Fig. 3. Secure operation of the grid imposes that, except for the RoCoF at COI, the RoCoF at the CP of cluster i to the main grid should also adhere to the specified 𝑅𝑜𝐶𝑜𝐹 𝑑𝑒𝑠 . More precisely, one can express the constraint on the local RoCoF at the CP of cluster i to the grid based on ( 3) and ( 6) as

∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 2𝐻 𝐷𝑉𝑃𝑃 𝑖 ≤ 𝑅𝑜𝐶𝑜𝐹 𝑑𝑒𝑠 → 𝐻 𝐷𝑉𝑃𝑃 𝑖 ≥ ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 2𝑅𝑜𝐶𝑜𝐹 𝑑𝑒𝑠 ( 10
)
The right-hand side of (10) specifies the minimum required inertia for cluster i to comply with the grid codes regarding RoCoF.

By calculating minimum 𝐻 𝑖 from ( 10), one can determine the contribution factor (CF) of cluster i from a dynamic perspective, denoted as 𝐶𝐹 𝑖 𝐷 , as

𝐶𝐹 𝑖 𝐷 = 𝐻 𝐷𝑉𝑃𝑃 𝑖 𝐻 𝐷𝑉𝑃𝑃 (11)
While ( 11) is used to calculate n-1 CFs related to the associated clusters, the n th CF D is calculated in a way that satisfies the overall 𝐻 𝐷𝑉𝑃𝑃 . To achieve this, one can write

𝐶𝐹 𝑛 𝐷 = 1 - ∑ 𝐻 𝐷𝑉𝑃𝑃 𝑖 𝑛-1 𝑖=1 𝐻 𝐷𝑉𝑃𝑃 ( 12
)
It is noteworthy that although equations ( 11) and ( 12) should produce identical results, measurement errors in RoCoF and simplifications in calculations can introduce some differences. This implies that the n th cluster is designated as the slack cluster in [START_REF] Yan | Data-driven load frequency control for stochastic power systems: A deep reinforcement learning method with continuous action search[END_REF] to compensate for the balancing power required to precisely meet 𝐻 𝐷𝑉𝑃𝑃 .

The physical interpretation of ( 11) and ( 12) holds true only when either the RoCoF at COI or the CP of clusters exceeds the desired value; otherwise, the CF D is regarded as zero.

Steady-state point of view: The steady-state frequency deviation following occurrence of the power imbalance ∆𝑃 𝐿 can be calculated by [START_REF] Bevrani | Robust power system frequency control[END_REF] ∆𝑓 𝑠𝑠 = ∆𝑃 𝐿 𝛽 𝐺 [START_REF] Sadeghi | Optimal bidding strategy of a virtual power plant in day-ahead energy and frequency regulation markets: A deep learning-based approach[END_REF] where

𝛽 𝐺 = 1 𝑅 𝐺 + 𝐷 𝐺 ( 14
)
In ( 13) and ( 14) ∆𝑓 𝑠𝑠 , 𝛽 𝐺 , 𝑅 𝐺 and 𝐷 𝐺 are steady-state frequency deviation, frequency response characteristic of the grid, droop characteristic of the gird, and load damping of the grid, respectively.

On the other hand, the steady-state frequency deviation in response to the contribution of DVPP in PFC can be explained using the equivalent frequency response characteristics i.e., 𝛽 𝐺 + 𝛽 𝐷𝑉𝑃𝑃 , based on (13), as

∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 =
∆𝑃 𝐿 𝛽 𝐺 + 𝛽 𝐷𝑉𝑃𝑃 [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] and

𝛽 𝐷𝑉𝑃𝑃 = 𝑟 𝑅 𝐷𝑉𝑃𝑃 + 𝑟𝐷 𝐷𝑉𝑃𝑃 ( 16
)
where 𝑟, ∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 and 𝛽 𝐷𝑉𝑃𝑃 are participation factors of DVPP for the under-study operating point, desired steady-state frequency deviation after provision of FRR, and aggregated actual frequency response characteristic of the DVPP, respectively.

While the aim is to calculate 𝛽 𝐷𝑉𝑃𝑃 in [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF], ∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 is also an unknown parameter. However, ∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 can be calculated from the SO's point of view where the DVPP can be represented by the FRR power ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 . More precisely, ∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 can be calculated as

∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 = ∆𝑃 𝐿 -∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 𝛽 𝐺 ( 17
)
By having ∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 from [START_REF] Ochoa | Control Systems for Low-Inertia Power Grids: A Survey on Virtual Power Plants[END_REF], one can calculate 𝛽 𝐷𝑉𝑃𝑃 for the overall DVPP, based on [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF], as

𝛽 𝐷𝑉𝑃𝑃 = ∆𝑃 𝐿 ∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 -𝛽 𝐺 ( 18
)
To calculate the contribution of each cluster of sources within the DVPP to realize 𝛽 𝐷𝑉𝑃𝑃 [START_REF] Yan | Frequency Control and Optimal Operation of Low-Inertia Power Systems with HVDC and Renewable Energy: A Review[END_REF], again consider the i th cluster in Fig. 3. Secure operation of the grid imposes that the steady-state frequency deviation at the CP of the clusters to the main grid should be the same as that of ∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 in [START_REF] Ochoa | Control Systems for Low-Inertia Power Grids: A Survey on Virtual Power Plants[END_REF]. More precisely, one can express the local steady-state frequency deviation in the CP of cluster i to the grid based on (3) and ( 13) as

∆𝑓 𝑠𝑠 𝑖 = ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 ∆𝑓 𝑠𝑠 𝑖 =∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 → 𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 = ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 ∆𝑓 𝑠𝑠 𝑑𝑒𝑠𝑖𝑟𝑒𝑑 ( 19
)
where ∆𝑓 𝑠𝑠 𝑖 and 𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 are steady-state frequency deviation at the CP of cluster i to the main grid and frequency response characteristic of cluster i, respectively. By calculating 𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 from [START_REF] Marinescu | POSYTYF concept and objectives[END_REF], one could calculate the CF of cluster i from a steady-state perspective, denoted as 𝐶𝐹 𝑖 𝑠𝑠 , as

𝐶𝐹 𝑖 𝑠𝑠 = 𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 𝛽 𝐷𝑉𝑃𝑃 (20)
Again, by considering the n th cluster as slack cluster, one can write

𝐶𝐹 𝑛 𝑠𝑠 = 1 - ∑ 𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 𝑛-1 𝑖=1
𝛽 𝐷𝑉𝑃𝑃 [START_REF] Lasseter | Grid-forming inverters: A critical asset for the power grid[END_REF] The final CF of each cluster, which satisfies both the dynamic and steady-state aspects, can be defined as

𝐶𝐹 𝑖 = max{𝐶𝐹 𝑖 𝑠𝑠 , 𝐶𝐹 𝑖 𝐷 } ( 22
)
With the CF for cluster i determined, the contribution of the cluster in realizing FRR power ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 , represented by ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 𝑖 , which aligns with the associated reference signal denoted as ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 𝑖

, is defined as

∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 𝑖 ≅ ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 𝑖 = 𝐶𝐹 𝑖 × ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 ( 23
)
2) AGC PFC stabilizes frequency at a value lower than the nominal one. Therefore, the SO needs to provide an auxiliary reference signal to the DVPP in order to restore frequency to the nominal value. To achieve this, the difference between ∆𝑃 𝐿 and ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 should be compensated by the ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 in Fig. 1. However, the primary concern is how to allocate the reference signal among the sources within the DVPP.

To disaggregate the reference signal, let's define the overall ACE as [START_REF]Balancing and Frequency Control Reference Document[END_REF] 𝐴𝐶𝐸 = 𝐴𝐼𝑃 -𝑆𝐼𝑃 -𝐵(𝑓 𝑎 -𝑓 𝑠 ) [START_REF] Adabi | Direct Participation of Dynamic Virtual Power Plants in Secondary Frequency Control[END_REF] where AIP, SIP, B, fa, and fs are actual interchanged power, scheduled interchanged power, frequency bias, actual frequency, and nominal frequency, respectively. As the steady-state frequency is common throughout an interconnection following provision of the FRR, one can write ACE for cluster i, based on [START_REF] Adabi | Direct Participation of Dynamic Virtual Power Plants in Secondary Frequency Control[END_REF], as

𝐴𝐶𝐸 𝑖 = 𝐴𝐼𝑃 𝑖 -𝑆𝐼𝑃 𝑖 -𝐵(𝑓 𝑖 -𝑓 𝑠 ) (25)
As PFC and AGC affect system dynamics via emulated turbine-governor behavior but using different inputs, the reference signal for cluster i, can be calculated using the same procedure as in [START_REF] Marinescu | Dynamic virtual power plant: A new concept for grid integration of renewable energy sources[END_REF], namely

∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖 = 𝐴𝐶𝐸 𝑖 𝐴𝐶𝐸 (∆𝑃 𝐿 -∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 ) ( 26
)
where ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖

denotes SFC reference signal for cluster i.

F. Real-time re-adjusting of the setpoints 1) PFC Now, let's assume that in response to changing the operating point or the size of power imbalance, the output power of sources changes. Eq. ( 27) illustrates that in order to obtain the updated reference signal, the new CF needs to be calculated. First, let to consider 𝐶𝐹 𝑖 = 𝐶𝐹 𝑖 𝑠𝑠 in [START_REF] Johnson | Synchronization of parallel single-phase inverters with virtual oscillator control[END_REF]. Considering the long-term nature of the frequency dynamics in frequency regulation studies, one could express the deviation of CFi based on the Taylor expansion form of ( 20 where ∆𝑝 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 refers to the pre-fault steady-state change in the output power of cluster i in response to a changing operating point.

According to [START_REF] Marinescu | POSYTYF concept and objectives[END_REF], the second term on the right-hand side of (28) could be re-written as [START_REF] Golpîra | Emulation of virtual inertia to accommodate higher penetration levels of distributed generation in power grids[END_REF] Following the same notation as in [START_REF] Golpîra | Emulation of virtual inertia to accommodate higher penetration levels of distributed generation in power grids[END_REF], the third term on the right-hand side of (28) could be represented by Substituting ( 31) into (30) using the definition of ( 15) gives

By replacing ( 29) and ( 32) into ( 28) and rearranging the terms, the final equation takes the form

Now, let to consider 𝐶𝐹 𝑖 = 𝐶𝐹 𝑖 𝐷 in [START_REF] Johnson | Synchronization of parallel single-phase inverters with virtual oscillator control[END_REF]. By following the same procedure as in [START_REF] Oureilidis | Ancillary services market design in distribution networks: Review and identification of barriers[END_REF], one can write

∆𝐶𝐹 𝑖 = 𝐻 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 ∆𝑅𝑜𝐶𝑜𝐹 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 (1 -∆(∆𝑃 𝐿 ))
∆𝑃 𝐿 [START_REF] Kundur | Power system stability and control[END_REF] The new contribution factor 𝐶𝐹 𝑖 𝑛𝑒𝑤 can be computed by replacing ∆𝐶𝐹 𝑖 , obtained from [START_REF] Oureilidis | Ancillary services market design in distribution networks: Review and identification of barriers[END_REF] or [START_REF] Kundur | Power system stability and control[END_REF] based on the selection made in [START_REF] Johnson | Synchronization of parallel single-phase inverters with virtual oscillator control[END_REF], into [START_REF] Golpîra | A center-of-gravity-based approach to estimate slow power and frequency variations[END_REF]. Indeed, [START_REF] Oureilidis | Ancillary services market design in distribution networks: Review and identification of barriers[END_REF] and [START_REF] Kundur | Power system stability and control[END_REF] demonstrate that the PFC setpoint of the sources can be readjusted by simply monitoring changes in the output power of the sources within the DVPP, the disturbance size, and the local RoCoFs.

2) AGC In another effort to fully achieve an online frequency regulation scheme, the SO is going to update the reference signal ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝐼 in [START_REF] Milano | Frequency divider[END_REF] in response to changing ACE. Following the same procedure as in ( 27)- [START_REF] Kundur | Power system stability and control[END_REF] The updated SFC reference signal for cluster i can be calculated by substituting [START_REF] Golpîra | A datadriven under frequency load shedding scheme in power systems[END_REF] in [START_REF] Zhang | Deloading power coordinated distribution method for frequency regulation by wind farms considering wind speed differences[END_REF].

G. Scalability

The realization of ( 33) through [START_REF] Golpîra | A datadriven under frequency load shedding scheme in power systems[END_REF] in an online application requires measuring two sets of measurements, including 1) local dynamics of 𝐴𝐶𝐸 𝑖 , 𝑅𝑜𝐶𝑜𝐹 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 and variations in the output power of resources (∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 ), and 2) a global dynamic 𝐴𝐶𝐸. As the dynamic grid model in frequency regulation studies is divided into coherent subsystems connected by a limited number of tie-lines, the global dynamic 𝐴𝐶𝐸 signal and its variations can be easily measured by monitoring the exchanged power through these tie-lines. Hence, it becomes evident that the proposed framework can alleviate the computational burden by compressing the required measurements in the tie-lines and a few local measurements, which demonstrate scalability of the method.

H. Technical Considerations and Limitations

Before validating the proposed control laws through time domain simulation (TDS) results, it is essential to address certain technical aspects concerning the involvement of IBRs at the distribution level in ancillary services. Given that the DVPP, as a new entity, actively engages in the ancillary service market, it focuses solely on constraints linked to its generation resources, overlooking technical constraints imposed by the distributed network. These include limitations such as line and transformer thermal thresholds, over/under voltages in both steady-state and dynamic scenarios, and the erroneous activation of protection mechanisms. This query could be addressed by incorporating the distribution system operator (DSO) into the contingency analysis conducted by the SO. To be more precise, the approved contract between the SO and DVPP must receive the DSO approval before implementation. The method presented in this paper only revolves around distributing regulation services among clusters of IBRs within DVPP. The exploration of the potential and challenges associated with DSO-SO coordination is discussed in [START_REF] Oureilidis | Ancillary services market design in distribution networks: Review and identification of barriers[END_REF].

III. SIMULATION AND RESULTS

Exploratory studies are conducted to assess the effectiveness of the proposed formulations using two test systems: 1) a two-area power system, and 2) a large-scale New York-New England (NYNE) test system.

A. Two-Area test system Figure 4 shows a modified single-line diagram of the twoarea test system, adapted from [START_REF] Kundur | Power system stability and control[END_REF]. The system consists of four SGs and a DVPP, realized using an array of PVs in the distribution level, connected to the main grid via transformer T1, and a wind farm in transmission level in area 2. The level of RESs penetration, denoted as r in Fig. 1, is established at 35%. Among these, 65% operate in grid-forming mode as dispatchable sources, while 35% operate in grid-following mode as non-dispatchable sources. The IBRs within the DVPP functioning in grid-forming mode are presumed to supply regulatory power, encompassing FRR and AGC. A regulation reserve for wind farms and solar panels is ensured by considering a deloading primary control. This provides active power margin in the maximum power point tracking (MPPT) algorithm, usually considered to be 10% of the nominal power of the generator [START_REF] Adabi | Direct Participation of Dynamic Virtual Power Plants in Secondary Frequency Control[END_REF][START_REF] Zhang | Deloading power coordinated distribution method for frequency regulation by wind farms considering wind speed differences[END_REF]. To evaluate the effectiveness of the proposed method, a 5% load disturbance representing a power surge of 58.4 MW is applied at bus 5 in area 1. It is assumed that according to the PFC market settlement, the DVPP is responsible for providing 20 MW FRR i.e., ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 = 20. The steady-state frequency deviation in response to the applied disturbance can be calculated based on [START_REF] Sadeghi | Optimal bidding strategy of a virtual power plant in day-ahead energy and frequency regulation markets: A deep learning-based approach[END_REF] as

∆𝑓 𝑠𝑠 = 0.05 0.5653 = 0.088 ( 37
)
The provision of the FRR by the DVPP imposes, as per ( 17), a steady-state frequency deviation of Now, let to calculate the CFs for the sources within the DVPP. For this purpose, the injected power of the PV arrays through T1 to the main grid versus measured frequency at T1, as shown in Fig. 5, is applied to the Prony analysis. Decomposing the resulting signal, as illustrated in Figure 5, using (4), leads to Furthermore, one can write for cluster 2 as:

𝐻 𝐷𝑉𝑃𝑃 2 = 0.1373 (41𝑏)
Taking into the measured RoCoF of 0.013 pu/sec in the CP of the PV arrays to the main grid and utilizing the power imbalance propagation paradigm (3)-( 5), the local size of the power imbalance can be calculated as ∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 = 0.41422 × 0.013 = 0.00538 (42)

By having the power imbalance size (42) and the desired steady-state frequency deviation (38), 𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 1 can be calculated based on [START_REF] Marinescu | POSYTYF concept and objectives[END_REF] as

𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 1 = 0.00538 0.0582 = 0.092 (43)
Finally, the CF of the PV arrays in PFC can be calculated using [START_REF] Marinescu | Dynamic virtual power plant: A new concept for grid integration of renewable energy sources[END_REF] as follows

𝐶𝐹 1 𝑠𝑠 = 0.092 0.294 = 0.31 (44)
From a dynamic perspective, the RoCoF at the COI can be calculated using [START_REF] Kumar | Decentralized reinforce-control of renewable dynamic virtual power plant enabling it to offer ancillary services: An attempt towards net-zero targets[END_REF], as The RoCoF obtained from (45) inherently meets the target value of 0.5 Hz/s without requiring any DVPP involvement, effectively eliminating the minimum DVPP contribution needed for inertial response. Moreover, the second derivative of voltage angle multiplied by 2π at the CP of clusters 1 and 2 also falls within the acceptable range, resulting in the CF of the clusters being zero as well. Therefore, one can write for the overall CF based on [START_REF] Johnson | Synchronization of parallel single-phase inverters with virtual oscillator control[END_REF] as

𝑅𝑜𝐶𝑜𝐹 = 0.05 2(7.394) × 60 = 0.2 [ 𝐻𝑧 𝑠 ] (45)
𝐶𝐹 1 = max{0.31, 0} = 0.31 (46)
Furthermore, one could write for the CF of the wind farm, based on (21), as

𝐶𝐹 2 = 1 -0.31 = 0.69 (47)
When calculating CF2 using (20), a discrepancy of 1.3% is observed compared to (47), primarily due to measurement errors.

The calculation for determining the PFC setpoint signal to be assigned to cluster 1 within the DVPP can be derived using [START_REF] Kontis | Provision of primary frequency response as ancillary service from active distribution networks to the transmission system[END_REF] in the following manner:

∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 1 = 0.31 × 20 = 6.2 (48)
Furthermore, one can write for cluster 2 as

∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 2 = 0.69 × 20 = 13.8 (49)
Now, let's compare the PFC results given in (48) and (49) with those obtained using [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF]. Table I provides the comparison of the contribution of clusters 1 and 2 in providing active power for PFC.

The findings demonstrate that while the approach of [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] achieves the contracted PFC with a negligible error, the proposed method precisely fulfills the specified 20 [MW] requirement. Of greater significance, both methods achieve frequency regulation with a single control action, indicating that the primary setpoint signals produced by the control laws are effectively meeting the desired dynamics. Now, let's explore a situation where the output power of cluster 1 is bounded to a maximum of 6 [MW]. The proposed method, once again, generates the outcomes presented in Table II using a single control action, while [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] requires three control actions to achieve convergence to the results. Reliance on only one control action, ensuring low computational burden, for the proposed method makes it attractive for practical applications. Nonetheless, the incorporation of the hard constraint into [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] poses a cut-off value to the control input, leading to discrepancies between the desired and actual frequencies after each control action. As a result, it required the control laws to generate new setpoints. When coordinating with SGs to provide PFC, these successive adjustments of the IBR setpoints lead to unnecessary governor mechanical movements, which reduces the governor's lifetime. On the other hand, the AGC system should compensate for the frequency deviation of (38). To accomplish this, one can express the reference signal for the PV arrays in Fig. 4, based on (26), as with the use of ( 26), a disparity of 0.9% is noted when compared to the value obtained from (51), mainly attributable to measurement errors.

Figure 6 compares the frequency responses (FRs) of the system with and without the contribution of the DVPP. The findings demonstrate that the PFC and AGC setpoint signals transmitted to clusters 1 and 2, as referenced in (48) via (51), are capable of arresting frequency decline and returning it to its nominal value. In this context, it is crucial to emphasize that the slight difference in frequency nadir between scenarios with only PFC and those incorporating both PFC and AGC can be attributed to the different number of IBRs required to operate in grid-forming mode to meet AGC requirements, as opposed to using PFC alone. More precisely, in response to the deloading of resources that contribute to SFC and to maintain the targeted penetration level of RESs in the system, the quantity of wind turbines has been increased, thereby accounting for this discrepancy. Now, let to consider a situation where the output power of the sources in clusters 1 and 2 changes in response to weather changes. The new PFC setpoint signal can be calculated using both [START_REF] Golpîra | A center-of-gravity-based approach to estimate slow power and frequency variations[END_REF] Recalculating the exact setpoint signal using the same procedure as in ( 37)-(49) results in the new setpoint of 4.54. This means that (52) yields a result with an error of approximately 0.44%, which is promising for online applications.

Another attempt is made to calculate the new setpoint signal in response to changes in both weather conditions and the applied fault. For his purpose, ∆𝐶𝐹 𝑖 can be written using [START_REF] Oureilidis | Ancillary services market design in distribution networks: Review and identification of barriers[END_REF] as The new PFC setpoint signal reported in (54) is then calculated by substituting (53) into [START_REF] Golpîra | A center-of-gravity-based approach to estimate slow power and frequency variations[END_REF], which is well fitted to the results obtained using the same procedure as in ( 37)-( 49), with an error of 0.31%. Moreover, it is expected that frequency trace for the case with the contribution of the DVPP experiences smaller frequency deviation at the nadir point compared to the trace without the DVPP. The variation in the nadir between both traces is attributed to the DVPP's contribution to the system's inertial response. To calculate this contribution, consider the derived equation in [START_REF] Golpîra | A datadriven under frequency load shedding scheme in power systems[END_REF], as

𝐻 = 𝑇 𝑛 ∆𝑃 𝐿 4∆𝑓 𝑛𝑎𝑑𝑖𝑟 ( 55
)
where 𝑇 𝑛 and ∆𝑓 𝑛𝑎𝑑𝑖𝑟 are frequency nadir time and variation of frequency nadir, respectively.

Referring to the traces associated with frequency response with and without the contribution of DVPP in Figure 6, and utilizing equation ( 55), one can express the following: The mathematically derived inertia deviation of (56) aligns closely with the calculated inertia constant for the DVPP in (41a) and (41b). This further substantiates the accuracy of the proposed methods in this paper. Now, let to assume that the applied disturbance in the system changes by 10 MW i.e., ∆𝑃 𝐿 = 48.4. The variation in the SFC reference signal for cluster 1 can be calculated using [START_REF] Golpîra | A datadriven under frequency load shedding scheme in power systems[END_REF] To enable precise compensation for the change in power imbalance within the AGC system, one can express the variation of the setpoint signal for the slack cluster as ) using [START_REF] Golpîra | A datadriven under frequency load shedding scheme in power systems[END_REF], a discrepancy of 1.2% is observed in comparison to the value obtained from (58).

∆𝐻 = 2.
B. NYNE test system In another attempt, the efficiency of the proposed method is investigated on a large-scale system. Figure 7 presents a single-line diagram of the NYNE test system, illustrating the generator buses and coherent areas. The system consists of five geographical regions, 68-Bus, 16-Machine, and a DVPP comprising 5 cluster of RESs. In this representation, generators G1 to G12 are equipped with a fourth-order type II power system stabilizer, tuned to provide sufficient damping [START_REF] Rogers | Power system oscillations[END_REF]. It is assumed that the SGs have governor in service. Additionally, the DVPP participate in both PFC and SFC services. The realization of the DVPP is as follows The RESs penetration level is set at 35%, with 70% functioning in grid-forming mode and 30% in grid-following mode. Of note that 15% of the nominal power of wind farms and 10% of the nominal power of solar panels functioning in grid-forming mode are considered as an active power margin for regulation service.

To evaluate the efficiency of the proposed method on the NYNE test system, four different scenarios are considered as ➢ Scenario 1: Tripping of G2 in area 1 without fault ➢ Scenario 2: Tripping of G10 in area 2 without fault ➢ Scenario 3: Switching on of a load at bus 7 ➢ Scenario 4: Shedding of 24.7 pu active load in bus 18 Firstly, the size of power imbalances in the CP of the clusters to the main grid for the aforementioned scenarios is calculated using power imbalance propagation paradigm. Table III reports the size of power imbalances, which are obtained by following the same procedures as those in (39)-(42).

To further evaluate the accuracy of the results in Table III, the quantities are compared with those obtained using local dynamics estimation methods [START_REF] Milano | Frequency divider[END_REF][START_REF] Golpîra | A center-of-gravity-based approach to estimate slow power and frequency variations[END_REF]. The maximum errors associated with each scenario are reported in Table IV.

The per unit PFC setpoints (∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 𝑖

), obtained through the identical methods outlined in (37) to (49), for the five cluster of sources within the DVPP, are reported in Table V. Moreover, the per unit reference signals transmitted by the AGC system to restore the frequency to its nominal value, obtained using (50) and (51), are reported in Table VI. depict the frequency responses of the center of inertia (COI) for the scenarios. The analysis of TDS results is employed to validate the design of control-responsive outputs allocation to IBRs. In this context, the allocation of reference signals among the resources within DVPP for participation in PFC is related to stability issues, ensuring that the system remains within the permissible range for RoCoF and frequency nadir. Based on Fig. 8 and the reported results in Table VII, it is apparent that, although the RoCoF for the trace without DVPP contribution deviates from the acceptable value of 0.5 Hz/sec, the proper distribution of FRR among the IBRs guarantees operation within the allowable RoCoF range. Furthermore, Table VII compares the values of RoCoF for the traces in Figs. 8-11 with those obtained using [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF], highlighting the high efficiency of the proposed method in enhancing the frequency behavior of the system. On the other hand, the TDS validates the allocation of resources to provide SFC, which is a regulatory duty rather than stability-related issue. In another analysis, the variation of PFC setpoint signals in response to changes in the weather and applied fault is calculated using [START_REF] Golpîra | A center-of-gravity-based approach to estimate slow power and frequency variations[END_REF] and either [START_REF]Balancing and Frequency Control Reference Document[END_REF] or [START_REF] Oureilidis | Ancillary services market design in distribution networks: Review and identification of barriers[END_REF], with the specific equation selection governed by [START_REF] Johnson | Synchronization of parallel single-phase inverters with virtual oscillator control[END_REF]. The variations in these parameters, represented as percentages of the primary values on which Table V The discrepancies between the values in Table IX and the exact results obtained using the same procedures as in ( 37)-(49) are also documented in Table X. These findings hold promise for the practical implementation of online applications for dispatching DVPP resources and facilitating frequency regulation. Finally, the effectiveness of ( 35) and [START_REF] Golpîra | A datadriven under frequency load shedding scheme in power systems[END_REF] 24)- [START_REF] Milano | Frequency divider[END_REF].

To further justify the scalability of the proposed approach, let's consider a scenario where the NYNE test system incorporates 5 DVPPs, each tasked with delivering ancillary services to specific areas. The SO has mandated these DVPPs to collectively provide a maximum FRR capacity of 123.5 MW, with individual allocations of 25 MW, 25 MW, 20 MW, 35 MW, and 18.5 MW for each DVPP, respectively. Now, imagine a situation in which the SO calls upon these DVPPs to deliver 85 MW of FRR in response to a load being switched on in area 1. Table XII compares the results of the proposed method with those obtained using [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] where a Lyapunov function is integrated into the RL-based method to optimize the allocation of PFC among resources. Although both methods can meet the system's FRR requirement, the proposed approach achieves setpoint adjustments in a single try, while [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] requires three iterative attempts.

Let's consider a scenario where, in response to real-time weather changes affecting the operation of wind farms operating in grid-following mode, the SO updated the command reference signal to 88 MW. Table XIII presents a comparison between the outcomes achieved through the suggested approach and those derived from [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF].

The outcomes presented in Table XIII demonstrate that minor fluctuations in the input data of the nonlinear RL-based controller from [START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] have a noteworthy impact on the reference signal for each resource involved in PFC service. When collaborating with SGs for PFC, these substantial variations in IBRs result in unnecessary mechanical movements of SG governors, thereby diminishing their operational lifespan. On the contrary, the approach proposed in this paper effectively addresses algorithmic stability by utilizing Taylor expansion around the operating point, ensuring moderate variations in the reference signals in response to small variations in input data.

IV. CONCLUSION

This research contributes to the development of an effective and practical approach for online frequency regulation considering the contribution of IBRs, offering valuable insights for power system operation and control. The study examines RESs simultaneously at both the distribution and transmission levels. Without loss of generality, the RESs are aggregated into a unified entity known as the DVPP, which is treated from the SO's perspective as a participant in the frequency regulation service. The proposed mathematically derived formulations aim to effectively allocate the FRR and AGC signals between the resources within the DVPP to restore the frequency to its nominal value. The results demonstrate that the proposed frameworks successfully regulate the frequency in real-time using PFC and SFC loops. Moreover, the study highlights the importance of designating one of the clusters as the slack cluster, which compensates for measurement errors and simplifications in the calculations.

While the TDS validates the allocation of resources to provide SFC, the internal dispatch of resources within DVPP should take into account the economics of resource dispatching. Furthermore, coordinating between DSOs and SOs is another area of interest for the successful integration of IBRs at the distribution level into regulatory services. These aspects should be further discussed in future works.

Fig. 1 .

 1 Fig. 1. Block diagram representation of frequency regulation service using SGs and DVPPs. The red blocks constitute the dynamics of the DVPP. Highorder non-linear models are used in simulations and simplified models are used for analysis; CF stands for contribution factor.where ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 , ∆𝑃 𝐷𝑉𝑃𝑃 𝑆𝐹𝐶 , ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶

Fig. 2 .

 2 Fig. 2. A general description of the proposed frequency regulation scheme.

Fig. 3 .

 3 Fig. 3. Power system including DVPP.

  ∆𝑃 𝐿 2(𝐻 𝐺 + 𝐻 𝐷𝑉𝑃𝑃 )≤ 𝑅𝑜𝐶𝑜𝐹 𝑑𝑒𝑠[START_REF] Dall'anese | Optimal regulation of virtual power plants[END_REF]

Fig. 4 .

 4 Fig. 4. Modified two-area test system.

Fig. 5 .

 5 Fig. 5. Injected power of cluster of interest to the main grid.

  By substituting the constitute components of the signal in (40) into (5), one could write the emulated inertia constant of the PV arrays, as seen from the gird, as 𝐻 𝐷𝑉𝑃𝑃 1 = 0.22 + 0.11 + 0.073 + 0.011 = 0.41422 (41𝑎)

Fig. 6 .

 6 Fig. 6. Comparison of FR of the two-area test system considering contribution of DVPP in frequency regulation.

  reference signal for cluster 2, designated as slack, can be expressed as

  10 -(-5.324) = -4.676 (58) After conducting the calculation of ∆(∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 2

Fig. 7 .

 7 Fig. 7. NYNE test system.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Comparison of FR of the NYNE test system with and without contribution of DVPP for scenario 1.

Fig. 10 .

 10 Fig. 10. Comparison of FR of the NYNE test system with and without contribution of DVPP for scenario 3.

Fig. 11 .

 11 Fig. 11. Comparison of FR of the NYNE test system with and without contribution of DVPP for scenario 4.

  The main problem of interest is to update ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 𝑖 in (23) based on a few measured quantities from the grid to meet the new PFC requirements, namely

	∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 𝑖,𝑛𝑒𝑤 = (𝐶𝐹 𝑖 + ∆𝐶𝐹 𝑖 ) ⏟	× ∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶	(27)
	𝐶𝐹 𝑖 𝑛𝑒𝑤

  ) as

	𝐶𝐹 𝑖 + ∆𝐶𝐹 𝑖 = 𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 𝛽 𝐷𝑉𝑃𝑃 + 𝜕𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 + 𝜕∆𝑃 𝐿 𝛽 𝐷𝑉𝑃𝑃 -1 𝜕𝛽 𝐷𝑉𝑃𝑃 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 𝜕∆𝑝 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 𝛽 𝐷𝑉𝑃𝑃 𝜕𝛽 𝐷𝑉𝑃𝑃 𝜕∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 ∆𝑝 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 𝑖 𝛽 𝐷𝑉𝑃𝑃 𝛽 𝐷𝑉𝑃𝑃 2 ∆(∆𝑃 𝐿 )	(28)

  , one could write

	∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖,𝑛𝑒𝑤 = ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖	+ ∆(∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖	)	(35)
	Applying Taylor expansion to (26) results in
	∆(∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖	) = = -	𝜕∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖 𝜕𝐴𝐶𝐸 𝑖 (∆𝑃 𝐿 -∆𝑃 𝐷𝑉𝑃𝑃 ∆(𝐴𝐶𝐸 𝑖 ) + 𝑃𝐹𝐶 ) ∆(𝐴𝐶𝐸 𝑖 ) 𝜕∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖 𝜕𝐴𝐶𝐸 𝐴𝐶𝐸 𝐴𝐶𝐸 𝑖 × (∆𝑃 𝐿 -∆𝑃 𝐷𝑉𝑃𝑃 𝑃𝐹𝐶 ) 𝐴𝐶𝐸 2 ∆(𝐴𝐶𝐸)	∆(𝐴𝐶𝐸) (36)

TABLE I COMPARISON

 I OF PFC FOR THE PROPOSED METHOD AND[START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF]

	Cluster	1	2
	Proposed method	6.20	13.80
	Ref. [15]	6.31	13.65

TABLE II COMPARISON

 II OF PFC FOR THE PROPOSED METHOD AND [15] CONSIDERING

	HARD CONSTARINT
	Cluster	1	2
	Proposed method	6	14
	Ref. [15]	6	13.98

  and[START_REF] Oureilidis | Ancillary services market design in distribution networks: Review and identification of barriers[END_REF], specifically

	∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 1.𝑛𝑒𝑤 = (0.31 +	1 0.294	-0.0014 0.0582 ) × 20 = 4.56 (52)

  • Wind farms in transmission level in bus 8, in area 1. • PV arrays in distribution level in bus 35, in area 2. • Wind farms in transmission level in bus 47, in area 3 • Wind farms in transmission level in bus 42, in area 4 • Wind farm in transmission level in bus 50, in area 5

  is based, are reported in Table VIII. Table IX reports the new PFC setpoint signals.

TABLE III SIZE

 III OF POWER IMBALANCES IN THE CP OF THE CLUSTERS TO THE MAIN GRID

	#	∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 1	∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 2	∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 3	∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 4	∆𝑃 𝐿 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 5
	1	0.00713	0.00321	0.00593	0.00538	0.00661
	2	0.00215	0.00391	0.00281	0.00407	0.00112
	3	0.00553	0.00641	0.00219	0.00201	0.00253
	4	0.00180	0.00012	0.00222	0.00103	0.00098

TABLE VI

 VI

		AGC REFERENCE SIGNALS IN PU
	∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖	Area 1	Area 2	Area 3	Area 4	Area 5
	Scenario 1	0.41	0.13	0.23	0.11	0.12
	Scenario 2	0.18	0.39	0.21	0.07	0.15
	Scenario 3	0.29	0.28	0.11	0.19	0.13
	Scenario 4	-0.22	-0.14	-0.33	-0.21	-0.10

  in updating the reference signals in response to changes in the ACE signals is demonstrated in Table XI. It is assumed that for each scenario in Table VIII, the size of the power imbalance decreases by 10%. The results reported in Table XI, illustrating the variation of the reference signals consistent with the findings in Table VIII, closely align with the precise results obtained by following the same procedures as in (

  ∆𝑃 𝐿 ) ∆𝑝 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 1 ∆𝑝 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 2 ∆𝑝 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 3 ∆𝑝 𝑐𝑙𝑢𝑠𝑡𝑒𝑟 4 ∆𝑝𝑐𝑙𝑢𝑠𝑡𝑒𝑟 5

				TABLE VII
	COMPARISON OF ROCOF RELATED TO COI FOR THE PROPOSED METHOD AND
				[15]
			Scenario	1	2	3	4
		Proposed method	0.0571 0.0413 0.0531 0.0163
			Ref. [15]	0.7619 0.0416 0.0531 0.1333
		Without DVPP	0.8315 0.0423 0.0534 0.2114
				TABLE VIII
			PARAMETERS VARIATIONS
	# ∆(1 -10%	5%	6%	-4%	3%		2%
	2	15%	-3%	8%	7%	6%		-7%
	3	15%	-4%	-2%	4%	10%		7%
	4	20%	-4%	9%	6%	-13%		11%
				TABLE IX
	NEW PFC SETPOINT SIGNALS IN RESPONSE TO CHANGING OF WEATHER
			CONDITION AND APPLIED FAULT
	∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 𝑖	Area 1	Area 2	Area 3	Area 4		Area 5
	1		0.23-0.04	0.11-0.07	0.20+0.15	0.08-0.03	0.38-0.01
	2		0.18+0.06	0.18-0.09	0.24-0.11	0.10-0.05	0.30+0.19
	3		0.06+0.17	0.45+0.06	0.23-0.10	0.15-0.09	0.11-0.04
	4		0.22+0.08	0.13-0.04	0.31-0.18	0.16+0.18	0.18-0.04

TABLE X ERROR

 X BETWEEN THE RESULTS OF

  TABLE IX AND THE EXACT VALUES ∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑃𝐹𝐶 𝑖

		Area 1	Area 2	Area 3	Area 4	Area 5
	1	1.8%	0.65%	~0.00%	0.43%	0.3%
	2	0.3%	0.42%	0.72%	0.81%	~0.00%
	3	0.87%	0.13%	~0.00%	~0.00%	0.37%
	4	0.12%	~0.00%	0.42%	~0.00%	~0.00%

TABLE XI VARIATION

 XI OF THE REFERENCE SIGNALS IN RESPONSE TO CHANGES IN THE

			ACE SIGNAL
	∆𝑃 𝐷𝑉𝑃𝑃-𝑟𝑒𝑓 𝑆𝐹𝐶 𝑖	Area 1	Area 2	Area 3	Area 4	Area 5
	Scenario 1	0.37	0.11	0.22	0.09	0.11
	Scenario 2	0.17	0.36	0.19	0.06	0.12
	Scenario 3	0.23	0.27	0.10	0.18	0.12
	Scenario 4	-0.26	-0.16	-0.34	-0.23	-0.11

TABLE XII COMPARISON

 XII OF FRR FOR THE PROPOSED METHOD AND[START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] FOR 85 MW

		LOAD DISTURBANCE
	DVPP	1	2	3	4	5
	Proposed method	22.8	15.2	13.5	15	18.5
	Ref. [15]	21.28	15.22	13	17	18.5

TABLE XIII COMPARISON

 XIII OF FRR FOR THE PROPOSED METHOD AND[START_REF] Cui | Reinforcement learning for optimal primary frequency control: A Lyapunov approach[END_REF] FOR 88 MW

		LOAD DISTURBANCE
	DVPP	1	2	3	4	5
	Proposed method	23.6	15.6	14.2	15.3	19.3
	Ref. [15]	9.71	24.61	7.41	27.77	18.5