A Novel Allocation Strategy Based on the Model Predictive Control of Primary Frequency Regulation Power for Multiple Distributed Energy Storage Aggregators

: As the amount of distributed energy storage (DES) in a power system continues to increase, it will not be long before there are multiple DES aggregators participating in frequency regulation, and the realization of their coordinated control is a critical topic of current research. This study focused on the primary frequency regulation (PFR) power allocation strategy among multiple DES aggregators participating in PFR. This study ﬁ rst calculated the PFR demand according to a system frequency response model of the power system with DESs. Next, a PFR power allocation model of DES aggregators was developed based on model predictive control. The objective of this model was to minimize the overall frequency regulation cost while satisfying all of the constraints of DESs. Finally, the distributed interior point method was used to solve the model rapidly. The correctness and e ﬀ ectiveness of the proposed model and algorithm were veri ﬁ ed on two uni ﬁ ed transmission and distribution systems with DES aggregators used to supply the PFR service. The results revealed that the proposed model could e ﬀ ectively allocate PFR power to the various types of energy storage, with the additional bene ﬁ ts of slowing down the shift in the state of charge for energy storage units and ensuring the continuity of energy storage participation in frequency regulation.


Introduction
With the increasing penetration of wind, photovoltaic, and other renewable energy sources into power grids, the equivalent inertia of power systems decreases gradually [1]. The current frequency regulation (FR) method based on conventional generators, however, cannot fully satisfy the FR demand of the power system. Therefore, investigating novel methods of FR is critical [2].
Due to its bi-directional regulation and fast response characteristics, energy storage has been increasingly used to provide FR service in power systems [3]. Generally, largescale energy storage can participate in FR by forming energy storage stations with hundreds of megawatt hours. Nevertheless, it is uneconomical to directly construct energy storage stations due to their expensive investment costs. Instead, with the rapid growth of grid-connected distributed energy storage (DES) in distribution networks, the FR ability of numerous DESs could be utilized effectively [4].
Numerous studies have been conducted to address the problems associated with the participation of DESs in FR. Ghazavi et al. proposed several frequency support strategies for DESs under low-grid-inertia conditions in [5]. In [6], a profit-maximizing strategy, including the FR phase and state of charge (SOC) recovery phase, was presented for the DES aggregator to participate in the primary frequency regulation (PFR) market. Considering the uncertain behavior of users, Sun et al. [7] proposed a dynamic DES control strategy based on the risk-averse multi-armed bandits method. In [8], a fast hierarchical control framework was designed for DESs to provide frequency and voltage support using the coordinated injection of active and reactive power. This framework addressed two hierarchical schemes, namely, physical and logical schemes. The physical scheme is designed to prioritize the power injection from DESs, while the logical scheme refers to the computation for the power commands of DESs concerning their operational constraints. Considering the FR requirements of the power system, an independent AGC control strategy based on the area control error signal was proposed in [9] for DESs to participate in secondary frequency regulation (SFR). The FR power allocation strategies employed in the above studies are based on the available regulation capacity of DES, which is characterized by a simple control mode, but does not take into account the differences in the regulation costs and operational constraints of different DESs. Model predictive control (MPC) is an optimal control technique that can solve optimization problems under multiple constraints; as a result, it is widely used in various aspects of power systems [10,11]. Aiming at the SFR jointly carried out by DESs and thermal power units, Wang et al. [12] proposed an optimal control method based on model predictive control (MPC). Considering the uncertainty of system operation, a two-layer MPC-based DES control strategy was developed in [13] to provide efficient control signals that respond to the SFR requirements of the power system. In the above studies, the management of DESs was unified by the power dispatching and control center [5,9,12,13] or an aggregator [6][7][8]. However, with the increase in the number of DESs, multiple DES aggregators may participate in FR, where the realization of their coordinated control is a critical topic of current research.
The control methods of power systems discussed in the emerging literature can be generally categorized into centralized and distributed categories [14]. Distributed control strategies are more suitable for the application of DES aggregators participating in PFR, as this type of approach can protect the private information of DESs. To date, distributed algorithms have been widely used in other fields of power grids, such as economic dispatch [15,16], optimal power flow [17,18], and reactive power optimization [19,20], while distributed methodologies were not utilized in FR problems until recently. In [21], the alternating direction multiplier method (ADMM) was adopted to solve the power allocation problem of joint FR with wind turbines and energy storage. To collectively provide ancillary FR service using distributed resources, a distributed control method based on consensus control protocol was designed in [22]. In [23], a distributed control strategy based on a finite-time consensus algorithm was proposed for the DES aggregator to provide an FR service. In [24], Xiong et al. developed a novel strategy for coordinating wind turbines and energy storage to achieve consensus control by employing a specified time consensus approach. In [25], the authors proposed a two-level synergetic unit commitment model to integrate FR into the unit commitment scheduling strategy and solved the model by utilizing the analytical target cascading algorithm. However, the aforementioned distributed algorithms typically require more iterations to converge, which may lead to the potential risks of action failure or frequency instability in response to the frequency command.
The distributed interior point method (DIPM) is an optimality condition decomposition method [26], whose essence is the distributed solution of the Karush-Kuhn-Tucker (KKT) optimality conditions of coupled constraints. Unlike the distributed algorithms adopted in [21][22][23][24][25], the DIPM relies on the second-order optimization method (Newton method) and hence enjoys superior convergence properties. Due to its strong convergence, the DIPM has been widely used in the multi-region coordinated optimization of power systems. In [27], the DIPM is applied to solve for a multi-region optimal reactive power flow. However, the proposed method requires a central coordinator to solve the correction equations of the coupled variables. Based on [27], Lu et al. [28] proposed a fully decentralized distributed interior point method without a central controller and applied it to the optimal power flow problem. Considering the risk management of a multi-regional interconnected energy system, Zhang et al. [29] constructed a risk-averse model of regional energy service providers and applied the DIPM to solve it. At present, the distributed interior point method is mainly applied to the optimal dispatch and has not been used to solve the frequency regulation problem.
In sum, the participation of DESs in primary frequency regulation (PFR) mainly has to contend with two problems. The first is that DESs in a power grid are currently managed by a dispatching center or an aggregator, and the power allocation problems that arise from multiple aggregators are not considered. The second is that in order to protect the privacy of aggregators, the power allocation model should be solved by a distributed algorithm, but the distributed algorithms that are currently applied to the FR problem require many iterations to converge. This study focuses on these problems and proposes a novel power allocation strategy for multiple DES aggregators participating in PFR and applies the DIPM to solve the model in a distributed manner.
The main contributions of this study are as follows: ( (2) The DIPM was applied to solve the PFR power allocation model in a distributed manner, which can ensure the speed and accuracy of the solution while simultaneously protecting the private information of each aggregator. Unlike the distributed algorithms adopted in [21][22][23][24][25], the DIPM has a second-order convergence speed.
The remainder of this paper is organized as follows. In Section 2, the calculation for the PFR demand considering DES aggregators is introduced. Then, the PFR power allocation model is described in Section 3. In Section 4, the distributed solution for the PFR power model is elaborated upon. The case studies and analytical results are presented in Section 5. Finally, the conclusion of this study is given.

Calculation of the PFR Demand
In Figure 1, the architecture of the PFR service provided by multiple DES aggregators is displayed. Under scenarios of power disturbances, such as generator failure and tie line disconnection, considerable frequency fluctuations in power grids could occur. Then, the power dispatching and control center could calculate the PFR demand of the entire system based on the power disturbance and send it to the controlled DES aggregators. With the help of the developed distributed optimization algorithm, each aggregator can determine its optimal power allocation plan and allocate the active power output increment reference to each DES to adjust its outputs, which ultimately contributes to ensuring the frequency stability of the entire system.  In this study, based on the conventional low-order frequency response model [30], DESs were incorporated to construct the system frequency response model participating in PFR, as shown in Figure 2. This model consists of three types of transfer functions: a reheat turbine unit, a DES, and a rotor motion equation of the equivalent synchronous generator to describe system frequency variation.
. . . For the reheat turbine unit, the most significant time constant is the turbine equivalent inertia time, identified as TR in Figure 2. This constant tends to dominate the response of the largest fraction of power output. Compared with TR, all of the other reheat turbine unit time constants are insignificant. Ignoring all of the other constants, the transfer function of the reheat turbine unit can be written as

MPC-based PFR power allocation
where Km is the mechanical power gain coefficient, FH is the ratio of total power generated by the high-pressure turbine, and R is the governor regulation coefficient. DES has the advantages of a rapid response, flexible controllability, and stable operation. A first-order inertia link is usually adopted to represent the frequency response of DES [31]. The transfer function of the DES model is where TB is the response time of the ith DES. The rotor motion equation used to describe the frequency variation in the equivalent synchronous generator of the power system is written as where D is the load-damping coefficient, H is the system equivalent inertia time, ΔPd is the power disturbance of the system, ΔPG refers to the turbine mechanical power increment, ΔPBi is the active power output increment of the ith DES, and Δf is the frequency deviation.
To achieve fast PFR, energy storage typically uses droop control to adjust the active power output according to the droop coefficient. However, the system characteristics should be considered when adjusting the droop coefficient. Therefore, the total droop coefficient of DESs should be optimized according to the system frequency response model.
According to the frequency response model with DESs participating in PFR, the frequency domain expression of the frequency deviation can be derived [32]: By applying the inverse Laplace transform, the time domain expression of the frequency deviation can be obtained as follows: Let the partial derivative of Equation (7) with respect to time t be 0; then, the time (tnadir) corresponding to the maximum frequency deviation can be obtained and the maximum frequency deviation can be deduced as follows: According to the requirements of the technical guidelines for the safety and stability control of power systems in China, the maximum frequency deviation was set to 0.5 Hz. Thus, the total droop coefficient of DESs can be calculated from Equation (8) After obtaining the PFR demand, the power dispatching and control center then sends it to each DES aggregator for further calculation.

Primary Frequency Regulation Power Allocation Model
As mentioned previously, this study developed a novel MPC-based power allocation model that incorporates both a predictive model and the receding horizon optimization model. In this section, the predictive model of aggregators is introduced first, followed by an elaboration of the receding horizon optimization model.

Predictive Model of Aggregators
DES has fast response characteristics, and its dynamic process can be expressed using the following first-order differential equation: where a ∈ A, and A is the set of DES aggregators; a Bi T is the response time of the ith DES in aggregator a; and The SOC of the DES will change during the charging and discharging processes. The transfer process of the SOC of the DES can be described as follows: , where , Therefore, the predictive model in Equations (10) and (11) for the ith DES in aggregator a can be expressed in the standard form as follows: , and matrix 1 1 Assuming that Na represents the number of DESs that are owned by aggregator a, according to Equations (12) and (13), the predictive model of aggregator a can be expressed as follows: Based on the continuous model (14) and (15), the discrete model with sampling time Ts can be obtained as follows:

Receding Horizon Optimization Model
(1) Objective The objective of the PFR power allocation problem is to minimize the FR cost of the DES participating in the PFR, which includes the active power deviation cost and the SOC offset cost: (2) Constraints The active power output increment of DES should satisfy the following constraint: denote the active power output increment at time k + n and its maximum value for the ith DES in aggregator a, respectively. The SOC of DES should satisfy its upper and lower limits as follows: where ， min refer to the active power output increment references of aggregator a and its ith DES at time k + n, respectively, and Δf (k + n|k) denotes the frequency deviation at time k + n. The sum of the active power output increment references of all aggregators should be equal to the PFR demand as follows:

Distributed Solution of PFR Power Allocation Model
For the convenience of discussion, the PFR power allocation model with multiple DES aggregators expressed by Equations (18)-(22) can be expressed as follows: (25) where x a = [x a(I) ; x a(C) ] represents the decision variables of aggregator a, x a(I) represents the internal variables of aggregator a, and x a(C) represents the coupling variables between aggregator a and other aggregators. Equation (24) consists of Equations (21) and (22), and Equation (25) consists of Equations (19) and (20). In order to protect the privacy of each DES aggregator, the problem shown in Equation (20) can be transformed into a distributed quadratic programming problem and solved in a distributed manner using the DIPM [28]. The solution process is described below.

Construction of the Unconstrained QP Model
Considering the DES aggregator a as an example, its PFR power allocation model can be rewritten as follows: (1) Construction of the Lagrange function By introducing slack variables m a and Lagrange multipliers (or dual variables) y a and z a to Equations (27) and (28), the Lagrange function can be obtained as follows: where m a ≥ 0, z a ≤ 0, r is the number of inequality constraints, and µ a denotes the barrier factor of aggregator a.
(2) Derivation for the Karush-Kuhn-Tucker (KKT) optimality conditions The necessary condition for the existence of a minimal value of the Lagrangian function (29) is that the partial derivatives of the Lagrangian function with respect to all variables are 0. This leads to the following system of KKT equations: Based on Equations (33) and (34), the following equation can be obtained: (3) Derivation of the correction equation The KKT optimality conditions (30)-(33) can be solved using Newton's method. The following correction equation can be obtained by retaining Δx a and Δy a and eliminating other variables: By merging the unconstrained QP models of all aggregators, the aggregated unconstrained QP model can finally be obtained as follows:

Approach for Solving the Unconstrained QP Model
The unconstrained QP model of aggregator a in Equation (39) can be written as a QP model, which is expressed in terms of the coupling variables' increments. Each aggregator first forms its QP model, and then obtains the aggregated unconstrained QP model by communicating with the other aggregators. Thereafter, each aggregator solves the aggregated unconstrained QP model independently to obtain the coupling variables' increments, which, in turn, updates the primal and dual variables. The specific process is outlined in what follows.
(1) The construction of QP models expressed in the coupling variables' increments For the unconstrained QP problem in Equation (39), the decision variables' increments are rearranged in the order of the internal variables' increments and coupling variables' increments, that is, Δx a = [Δx a(I) ; Δx a(C) ]. Then, corresponding adjustments are made to the dual variables' increments and the coefficient matrices, and the following equation can be obtained: where the coupling variables' increments are regarded as parameters. Therefore, Equation (41) can be written as follows: According to Equation (42), the internal variables' increments and dual variables' increments can be expressed in terms of the coupling variables' increments as follows: Substituting Equations (43) and (44) into Equation (41), the QP models expressed in the coupling variables' increments can be expressed as follows: where a 1 C denotes the coefficient matrix corresponding to the quadratic terms, a 2 C denotes the coefficient vector corresponding to linear terms, and a c is a constant.
(2) Solving for the coupling variables' increments After building the QP model (45)  Thus, the computational steps of the DIPM can be summarized as follows: Step 0: Give k = 0; set tolerance ε = 10 −4 ; for aggregator a, initialize its primal variables x a (k) and m a (k) and dual variables y a (k) and z a (k).
Step 1: Calculate complementary gap Gap a (k) for each aggregator and obtain the maximum complementary gap Gap max (k) by communicating with other aggregators.
If Gap max (k) < ε, then output the optimal solution and stop.

Case Studies and Analysis
Numerical experiments were conducted with the Dell Precision 3660 Tower with Intel ® Core TM i7-12700K CPU@3.60 GHz and 32 GB RAM. The algorithm was coded using MATLAB R2022a software.

Parameter Setup
To verify the effectiveness of the proposed method, a unified transmission and distribution system with two DES aggregators to provide PFR service was utilized, as shown in Figure 3. The parameters of the frequency response model of the system [34] are listed in Table 1. The parameters of DESs are listed in Table 2. The expected reference values of the SOC of DESs were set as 0.5. The sampling time was 0.05 s, the control cycle was 0.2 s, and the prediction horizon was 5 s. The convergence accuracy of the DIPM was set as ε = 10 −4 . Figure 3. Connection diagram of a unified transmission and distribution system with two distributed energy system aggregators. To illustrate the superiority of the proposed PFR power allocation strategy, it was compared with the traditional allocation strategy (i.e., the strategy based on available regulation capacity used in [5][6][7][8][9]) in both case 1 and case 2.

Result Analysis
To fully reflect the differences in the PFR power allocation strategy of various types of DESs, five representative DESs, that is, DES2, DES4, DES5, DES7, and DES8, were selected for detailed analysis.

Analysis of Case 1
The change in the system frequency after the system load increased by 45 MW abruptly is shown in Figure 4. After the 45 MW load step rise, when DESs do not participate in FR, the lowest system frequency was 49.26 Hz, and the maximum frequency deviation reached 0.74 Hz, which does not satisfy the frequency requirements for the safe and stable operation of the power system. When the DESs participated in FR, the maximum frequency deviation of the system was reduced to 0.5 Hz, and the steady-state frequency increased from 49.63 to 49.70 Hz.   Figure 5 shows the PFR power response of the aggregators under different allocation strategies in case 1. Aggregator 2 took on a heavier FR task than aggregator 1 with the proposed allocation strategy because its DESs had a lower FR cost coefficient. Adopting the traditional allocation strategy resulted in the same phenomenon, but this was because the available regulation capacity of aggregator 2 was larger than that of aggregator 1, meaning that its output was also larger than that of aggregator 1. Comparing the two allocation strategies, we could determine that the output of aggregator 2 with the proposed allocation strategy was larger than that of aggregator 2 with the traditional allocation strategy, while the opposite was true for aggregator 1. This was because the proposed strategy took into account the FR cost coefficients of the DESs; because the cost coefficient of aggregator 2 was lower, its output power was larger than that of aggregator 2 using the traditional allocation strategy. Similarly, aggregator 1 had a higher cost coefficient, and thus, its output was smaller than that of aggregator 1 using the traditional allocation strategy.  Table 3 shows the FR costs when using the two allocation strategies in case 1. Compared with the traditional strategy, the proposed allocation strategy took the FR cost into consideration. Therefore, when the proposed allocation strategy was adopted to allocate the FR power, aggregator 1, which had a higher FR cost coefficient, was allocated less FR power, resulting in its FR cost being lower than that when the traditional allocation strategy was used. Because the FR cost coefficient of aggregator 2 was lower, it took on more FR tasks, which made its FR cost higher. By comparing the total FR cost with the two allocation strategies, it could be determined that the proposed allocation strategy was more economical than the traditional strategy.  Figure 6 illustrates the outputs of DESs in case 1. Figure 6 reveals that DES8 had the largest output in the FR process since its active power deviation cost coefficient was the smallest, whereas DES2 had the smallest output due to it having the highest power deviation cost coefficient. The technical parameters of DES4, DES5, and DES7 were identical, but the initial SOC differed considerably. Among the three DESs, the output of DES4 was the smallest in the FR process because its initial SOC was less than the expected reference value of the SOC, whereas DES7 maintained the largest output because its initial SOC was larger than the expected reference value. The output of DES5 was between those of DES4 and DES7 because its initial SOC was the same as the expected reference value of the SOC.   Figure 7a reveals that the change in the SOC of DES8 was more evident than that of DES2 since DES8 exhibited the largest capacity and the largest output, while DES2 had the smallest capacity and the smallest output, which are presented in Table 2 and Figure 6. Although DES4, DES5, and DES7 had the same rated capacity, the change in their SOCs differed, as shown in Figure 7a-c. Among these three DESs with the same rated capacity, the larger the output was, the more the SOC changed.

Analysis of Case 2
After a sudden reduction in the system load by 40 MW, the change in the system frequency is illustrated in Figure 8.   Figure 8 reveals that after the load step-down disturbance occurred when the DESs participated in FR, the highest system frequency was 50.65 Hz, and the maximum frequency deviation was 0.65 Hz, which did not satisfy the frequency requirements for the safe and stable operation of the system. When the DESs participate in FR, the maximum frequency deviation of the system was reduced to 0.5 Hz, and the steady-state frequency was reduced from 50.33 to 50.28 Hz. Figure 9 shows the PFR power response of the aggregators under different allocation strategies in case 2. Similar to case 1, the DESs' active power deviation cost coefficient of aggregator 2 was lower than that of aggregator 1. Therefore, aggregator 2 had a larger output under the proposed allocation strategy. The output of aggregator 2 was also larger than that of aggregator 1 when adopting the traditional allocation strategy because aggregator 2 had a larger available regulation capacity. Compared with the two allocation strategies, it can be seen that the output of aggregator 2 with the proposed allocation strategy was larger than that when using the traditional allocation strategy, whereas the opposite was true for aggregator 1. The reason for this was that the proposed strategy considers the FR cost coefficient of the DESs and the cost coefficient of aggregator 2 was lower, and thus, its output power was larger than that of aggregator 2 using the traditional allocation strategy. Similarly, because aggregator 1 had a larger cost coefficient, its output was lower when using the proposed strategy than when using the traditional allocation technique.  Table 4 shows the FR costs under the two allocation strategies in case 2. Similar to case 1, because the proposed allocation method considers FR cost, for aggregator 1, which had a small FR cost coefficient, the FR cost using the proposed strategy was smaller than that when using the traditional strategy. Meanwhile, the FR cost of aggregator 2 under the proposed strategy was higher than that of the traditional strategy because aggregator 2 had a smaller FR cost coefficient and bore more FR tasks. Compared with the total FR cost, we could find that the proposed allocation strategy was more economical than the traditional allocation strategy.  Figure 10 shows the outputs of the DESs in case 2. Similar to case 1, DES8 undertook most of the FR tasks because of its smallest active power deviation cost coefficient, whereas DES2 had the largest active power deviation cost coefficient. Therefore, the charging power of DES2 was the smallest. The charging power of DES5 was larger than that of DES7 but smaller than that of DES4, as the DESs were in the charging state in case 2. Thus, under a sudden reduction of power load, for DESs with the same rated capacity, the smaller the initial SOC was, the larger the charging power was.  Figure 11 shows the change in SOC of DESs in case 2. Figure 11a reveals that although DES8 had the largest capacity, it had the largest charging power during the FR process, which resulted in its SOC change being more evident than that of DES2, which had a small capacity and less charging power. For DES4, DES5, and DES7, although their rated capacities were identical, the change in their SOCs differed, as shown in Figure 7a-c. For the three DESs with the same rated capacity, the larger the charging power was, the more the SOC changed.

Sensitivity of the PFR Power Allocation Model
The prediction horizon is an important parameter of MPC, and it is used to control how far into the future the MPC predicts the system response when evaluating the optimality of the control objective. A longer prediction horizon generally increases the control performance but also results in increasing the computational complexity. In this subsection, we present the analysis results of the sensitivity of the performance of the proposed model in case 1 with different prediction horizons, as shown in Table 5. As can be seen from Table 5, the optimization objective of the proposed model (i.e., the total cost of frequency regulation) decreased as the prediction horizon increased and the solution time increased. The optimization objective of the model decreased at an accelerated rate when the prediction horizon varied from 2 s to 5 s, but the decreasing speed slowed down considerably from 5 s to 10 s; meanwhile, the computation time consistently increased. Considering the optimization objective and the solution time, 5 s was the most appropriate prediction horizon for the model proposed in this study.

Performance of the DIPM
To demonstrate the effectiveness and superiority of the solution algorithm, the DIPM and other distributed optimization algorithms, such as ADMM and the synchronous alternating direction multiplier method (S-ADMM), were compared. Quadprog, which is a QP solver in MATLAB, was used as a benchmark for the comparison. Details on the computational performance in case 1 are presented in Table 6. The parameters of ADMM were set to the following: ρ = 0.05, λ0 = 1, and ε = 10 −4 . The parameters of S-ADMM were the same as that of ADMM. Table 6 reveals that the DIPM was consistent with Quadprog in terms of solution accuracy. From the perspective of computational efficiency, compared with ADMM, S-ADMM had a shorter solution time, but also had more iterations, which resulted in higher requirements for communication reliability; meanwhile, the DIPM was optimal in terms of iterations and solving time. The results revealed that the DIPM could guarantee a good solution speed while ensuring high accuracy. To demonstrate the scalability of the proposed allocation strategy, an extended system (i.e., a unified transmission and distribution system with four DES aggregators to provide PFR service) was established, as shown in Figure 12. The parameters of the frequency response model of the system can be found in [35], and the parameters of the DESs are listed in Table A1. To demonstrate the advantages of the proposed PFR power allocation strategy, a case in which the system load increased by 150 MW at t = 10 s was designed. According to Equation (9), the total droop coefficient KB of the DESs in this case was found to be −116.6. Table 7 shows the FR costs with the two allocation strategies applied in the extended system. It can be seen from the table that the power allocation strategy proposed in this paper was more economical than the traditional allocation strategy when applied to larger power systems or when there was an increased number of DES aggregators.  Table 8 shows the details of the computational performance of the four algorithms applied in the extended system. The parameters of the four algorithms were the same as those in Section 5.2.4. Comparing the results in Tables 6 and 8, it can be seen that when applied to a larger power system or when there was an increased number of DES aggregators, the iterations and solution time of ADMM and S-ADMM increased rapidly and their computational efficiencies could not meet the FR requirements of the power grid. Meanwhile, the DIPM only required nine iterations to converge and had a solution time of 0.1421s, which was much smaller than that of ADMM and S-ADMM. The computational performance in the extended system fully demonstrated the superiority of the DIPM for distributed solutions.

Conclusions
To solve the problem of insufficient FR capacity of power systems resulting from the high penetration of various sources of renewable energy, such as wind power and photovoltaic power, an MPC-based PFR power allocation model for DES aggregators was developed. Based on the simulation results, the following conclusions could be drawn: (1) The proposed PFR power allocation strategy can fully consider the differences in the power, capacity, SOC, and FR cost of various types of DESs and intelligently allocate the FR tasks. Furthermore, the proposed allocation strategy can slow down the SOC offset in the FR process to ensure the continuity of DES participating in FR. (2) By applying the DIPM technique for distributed solving, the information exchanged by each aggregator does not involve any private data or the variable information of a single DES, which can ensure the privacy of each aggregator. In addition, compared with ADMM and S-ADMM, the DIPM was superior in terms of both computational efficiency and regulation performance.
However, the proposed PFR power allocation strategy requires predicting the behavior of the system, and the implementation cost of this was relatively high. Although the DIPM was used to transform the model into an unconstrained quadratic programming problem for distributed solutions, its calculation requirements were still larger than those of traditional allocation strategies [5][6][7][8][9]. Moreover, data communication is an essential part of distributed algorithms in real-world scenarios. Future work will focus on the influence of partial information loss on the convergence of the algorithm in the data communication process and on finding a suitable correction strategy to make the algorithm more robust. Acknowledgments: We thank LetPub (www.letpub.com (accessed on 17 August 2023)) for its linguistic assistance during the preparation of this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.