Optimal Primary Frequency Reserve Provision by an Aggregator Considering Nonlinear Unit Dynamics

Power system operators are increasingly procuring frequency regulation reserves from distributed energy resources via aggregators. Provision of the frequency regulation services by the aggregators can be studied from different aspects and by using different models, ranging from very-short-run to long-term scheduling models. By developing a methodology for ensuring the operational model's compliance with the dynamic constraints, this article bridges the gap between economic studies of the aggregator's frequency containment reserve market participation on the one hand and the distributed energy resources' dynamic response on the other. The article presents a mixed-integer linear program for the aggregator's portfolio management that incorporates constraints of individual units' fast dynamics, and validates them via dynamic simulations. The results demonstrate the importance of including reserve activation constraints when modelling storage technologies and settling time constraints when modelling conventional units.

Abstract-Power system operators are increasingly procuring frequency regulation reserves from distributed energy resources via aggregators.Provision of the frequency regulation services by the aggregators can be studied from different aspects and by using different models, ranging from very-short-run to long-term scheduling models.By developing a methodology for ensuring the operational model's compliance with the dynamic constraints, this article bridges the gap between economic studies of the aggregator's frequency containment reserve market participation on the one hand and the distributed energy resources' dynamic response on the other.The article presents a mixed-integer linear program for the aggregator's portfolio management that incorporates constraints of individual units' fast dynamics, and validates them via dynamic simulations.The results demonstrate the importance of including reserve activation constraints when modelling storage technologies and settling time constraints when modelling conventional units.Index Terms-Aggregator, power system dynamics, primary frequency control, primary frequency reserve, frequency containment reserve, reserve market, day-ahead market.M. Miletić, H. Pandžić, and I. Kuzle are with the University of Zagreb Faculty of Electrical Engineering and Computing, 10000 Zagreb, Croatia (e-mail: marija.miletic@fer.hr;hrvoje.pandzic@ieee.org;igor.kuzle@fer.hr).
M. Krpan is with the University of Zagreb Faculty of Electrical Engineering and Computing, 10000 Zagreb, Croatia, and also with the Hitachi Energy Sweden, AB SE-771 80, Sweden (e-mail: matej.krpan@fer.hr).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TPWRS.2023

I. INTRODUCTION AND LITERATURE REVIEW
A GGREGATORS are intermediaries that introduce dis- tributed energy resources (DERs) to the wholesale electricity and reserves markets, i.e. primary, secondary and tertiary frequency response.We model primary frequency reserve which is in the European market setting known as the frequency containment reserve (FCR).Aggregators' participation in FCR markets has recently gained popularity both in the research modelling and the real-world applications.One of the first aggregators to participate in the FCR market managed to successfully prove its ability to use a fleet of large industrial loads to take part in the frequency control [1].However, most European countries are still creating or adjusting their market rules and prequalification procedures to allow the FCR market participation for aggregators as independent entities.Namely, an aggregator must pass a prequalification process as a single entity and prove that it can manage its aggregated DERs in a way that ensures a correct response to frequency deviations.For example, the Dutch transmission system operator (TSO) Tennet requires that aggregators prove their ability to successfully respond with full reserve activation to a specific set of disturbances within the prescribed time period (30 s or 2.5 minutes, depending on the test) and retain the power at the required level for 15 minutes after an activation [2].A detailed overview of market rules and technical requirements for FCR providers in Europe is presented in [3].A pre-qualified aggregator can bid up to its pre-qualified capacity in the market and is free to schedule the units accordingly.A question then arises: how much freedom does an aggregator have when providing an FCR response with heterogeneous units considering the TSO's FCR rules?
Analysing the literature on aggregators, two distinct approaches to modelling can be identified.The first one is focused on the optimal control algorithms for aggregators providing frequency regulation.Authors in [1] provided a control algorithm for industrial loads performing FCR.The presented algorithm is implemented on a large scale in France and the aggregator received a prequalification from the French TSO (RTE).An algorithm for an aggregator of thermostatically-controlled loads participating in the Nordic FCR-N market is described in [4].The authors developed scheduling and control algorithms and used a recurrent neural network to consolidate a large variety of thermal loads and enable the aggregator to issue control signals.Zhu and Zhang [5] developed an algorithm for coordination of batteries providing FCR with an emphasis on the state-of-energy (SOE) recovery.These papers are primarily focused on developing adequate real-time control systems for providing the FCR service, while the economic optimisation was secondary or neglected.
The second approach to the aggregator operation is the economic analysis, where authors try to find optimal bidding strategies for market-participating aggregators, considering various techno-economic constraints.An example of such approach is [6], where the authors performed a portfolio optimisation of an aggregator of residential heat pumps considering the devices' availability and reliability.In [7], the batteries providing FCR in the Nordic power system were analysed, while in [5] the optimal SOE range for FCR-providing batteries was determined.
Aggregators of electric vehicles are another segment in the literature, where their profitability is investigated considering the electric vehicles' arrival and departure times, and the batteries' SOE, see e.g.[8], [9].As a rule, the articles concerned with economic analysis feature units with constant droop gain settings over time and do not investigate the degrees of freedom the aggregator has when scheduling its units.
In general sense, the degrees of freedom are any variables that can be optimised by the aggregator, such as energy market and reserve market quantities and price bids, or technical parameters, while also considering the dimensions of these variables.In our case, we optimise the DA and FCR capacity bids, the latter through the units' droop setting.Our DA and FCR capacity bids are one-dimensional, while the droop settings of our units are four-dimensional, optimised with respect to each unit, time, DA and FCR market scenarios as well as FCR activation scenarios.
Recently, studies that combine optimisation with dynamic analysis have appeared, aiming to bridge the gap between them [10], [11], [12], [13], [14], [15].These studies focused on optimising frequency control services, but they do not comprehensively consider relevant nonlinearities: insensitivity, saturation, deadzone and rate limits.These nonlinearities need to be adequately taken into account according to the grid code requirements [2], [16] as they impact the response of FCR-providing units.
None of the nonlinearities were considered in [10], while [11] only considers saturation.Authors in [12], [14] consider saturation, but neglect insensitivity and deadband, and approximate the turbine governor response as linear (constant ramp rate), which is not always accurate as will be discussed later, especially on the time scale of the primary frequency control.Badesa et al. [13] consider only fast units described by first order dynamics with saturation, but neglect rate limits, deadzone and insensitivity.Guzman et al. [15] consider rate limits and saturation of fast energy storage providing secondary frequency control, but neglect combined effects of deadzone and insensitivity.
Furthermore, [10], [11], [12], [13], [14], [15] assume fixed ramp-rate of units, which might introduce errors, as elaborated in Section II-B.From an aggregator's perspective, its decisionmaking framework must be as accurate as possible in order for its service delivery to be grid compliant.
Finally, [10], [11], [12], [13], [14], [15] all propose frameworks from a system operator's perspective for determining frequency quality criteria and unit dispatching.The framework we propose is from an aggregator's perspective, i.e. a market player who does not have an insight into the grid topology, but wants to position itself in the FCR and DA markets while satisfying the imposed grid code requirements for FCR provision.To the best of our knowledge, there are no papers that discuss this approach.Furthermore, the novelty is that we tackle the issue of a nonlinear unit response by introducing the settling time constraints derived from dynamic models into the optimisation framework which is also something we have not encountered in the literature.
Considering the shortcomings of the above works, our research introduces a comprehensive mixed-integer linear programming (MILP) model of an aggregator with a diverse Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
portfolio of units participating in the DA and the FCR markets, taking into account dynamic responses of the aggregated units, as well as FCR activation.FCR activation is usually disregarded in the models, with an assumption that it is symmetrically activated and therefore positive and negative activations zero out in terms of the energy.This is an acceptable assumption when modelling generators only, but when any type of storage system with non-ideal efficiency is included in the model, reserve activation becomes nontrivial.We derive MILP constraints based on settling times which reflect specific dynamics of three different types of units.The MILP model decides the aggregator's bids in the FCR and the DA markets and distributes the FCR response among units by deciding on their droop gain settings.We use the model to answer the question how can the aggregators of heterogeneous DERs, using MILP, build and schedule their portfolios in a way that satisfies the technical requirements of the TSOs.
The original scientific contribution of this article is threefold: r A methodology for ensuring grid-code compliance of the mathematical programming model of an aggregator using a non-linear dynamic model of the aggregator's portfolio.
r A stochastic MILP model for an aggregator bidding in the FCR and the DA markets considering the reserve activation and remuneration, the droop gain limits, the rate limits as well as the analytically derived settling times of a diverse portfolio of units.
r A verification of the MILP model for grid-code compliance during the FCR activation using a non-linear dynamic model that considers the inherent controller insensitivity, the programmed deadzone, the rate limits and the saturation.The article is organised as follows.Section II describes the dynamic and the MILP models, Section III provides an overview of the data used in the study and describes the considered cases, Section IV presents the results, while Section V concludes the article.

II. METHODOLOGY
Our proposed methodology is summarised in Fig. 1.In the first step, the MILP is initialised with the input data and a basic set of constraints.The dynamic model represents the same set of units as the MILP, but in more details, and is therefore initialised with the same input data.In the second step, the MILP is solved and the results are fed as inputs to the dynamic model.The dynamic model receives as inputs the MILP-generated DA market schedule through variable p DA and the droop gain settings for the FCR provision through variable k.The next step checks the grid code compliance of the MILP-generated schedule by running the dynamic model and verifying the following conditions [2]: 1) Ability to activate the full volume of FCR within 30 seconds for positive and negative step changes of frequency; 2) Ability to maintain the activated FCR for at least 15 minutes; 3) Ability to activate the corresponding share of FCR for successive step-wise frequency changes and maintain the corresponding power output for at least 15 minutes in between the step changes; 4) Ability to achieve a linear activation of the full FCR volume within 2.5 minutes for an up-ramp and down-ramp of frequency deviation signal; 5) That the frequency measurements and power output during the above mentioned tests meet any other requirements regarding accuracy, sampling, combined effects of insensitivity and deadzone, etc.If the model is shown to be non-compliant with any of these conditions, the MILP constraints are adjusted.The process is repeated until the MILP model is compliant with all the conditions.We add new constraints to make sure the dynamic limits are included in the MILP, which may reduce the feasibility set, thus providing worse or the same objective function value.However, the feasibilty of the MILP model and covergence of our methodology is preserved if the dynamically-constrained model is feasible.Through this process, using the dynamic model described in Section II-C, we arrive to the MILP model described in Section II-D.
For the sake of completeness, we use the definitions of deadzone and insensitivity from the European network code on requirements for grid connection of generators [17], where deadzone (i.e.deadband) is "an interval used intentionally to make the frequency control unresponsive", while insensitivity is "the inherent feature of the control system specified as the minimum magnitude of change in the frequency or input signal that results in a change of output power or output signal".

A. Problem Description and Assumptions
The aggregator controls a set of DERs and bids on their behalf in the FCR and DA markets.Considering the DA and FCR market participation decisions, the aggregator signals the DERs to change their droop gain set-points (i.e.droop gains) so they deliver just the right amount of FCR, while satisfying the full activation time requirement.We build a model with the following assumptions: r The aggregator portfolio consists of four types of units: hydro turbine-generators (HTG), steam turbine-generators (STG), battery storage systems (BSS) and supercapacitor storage systems (SCS).We chose these four types of units since they represent a gamut of dynamic characteristics such as stored energy, response time and transient behavior.Furthermore, these units are deterministic and dispatchable compared to wind and PV which require special deloading schemes for ensuring reserve as well as power forecasting which is beyond the scope of this article.Since wind and PV are converter-interfaced, they can be represented by models similar to the ones we use for batteries and supercapacitors [18], [19].
r All units are under 10 MW which, under most European grid-codes, allows them to be connected to the distribution system.We assume that the DERs are distributed in a wide area and that their total installed capacity is several orders of magnitude lower than the size of the power system.This means that the impact of the aggregator's portfolio on the grid frequency dynamics is small, but it still must satisfy the imposed frequency response requirements per the gridcode [2], [16], [17].
r We neglect the impact of aggregator's portfolio on the grid power flows.This is justified because there are no technical constraints regarding grid power flow in the grid code requirements for FCR provision [2], [16], [17] and an aggregator has no insight into the real-time system topology and states beyond its units' point of common coupling.Concerns regarding power flow lay entirely on the system operator.
r DA and FCR market prices are uncertain and included in the model in the form of representative days.An additional source of uncertainty is the FCR activation, which is modelled as a different set of scenarios for each representative day.
r All units participate in the FCR market and all, except the supercapacitors, participate in the DA market as well.The aggregator is a price-taker in both markets.
r Energy storage units (BSS, SCS) must recover their SOE within 2 hours of FCR activation.This is in line with Regulation 2017/1485 [16].The SOE recovery is usually performed in the intraday market, but here, for simplicity, we assume that the recovery is performed the next hour after the FCR activation at the DA price.
r We do not consider the aggregator's remuneration policy towards the DER units, but take into account the units' operating costs.

B. Why Constant Ramp Rate Assumption is Inadequate in FCR Scheduling Problems
An aggregator represents a group of FCR-providing units, and the output power response of the FCR-providing units generally exhibits a nonlinear behaviour, which is a consequence of the underlying physical properties and control system design.This is contrary to the common assumption of constant ramp rates which result in linearly varying power output over time [11], [12], [13], [14].Although such assumption simplifies the optimisation problem formulation, in reality it may lead to erroneous conclusions because FCR-providing units are not behaving as assumed.
To illustrate the consequences of constant ramp rate assumption, the response of a laboratory-scale hydro turbine-generator at the Smart Grid Laboratory of the University of Zagreb Faculty of Electrical Engineering and Computing was recorded for two step changes in the power reference: from 1 to 6 kW and from 1 to 11 kW (ΔP 1 = 5 kW and ΔP 2 = 10 kW).The recorded measurements are provided in Fig. 2. The settling time for ΔP 1 is approximately 30 s (marked with the orange diamond), while the settling time for ΔP 2 is approximately 35 s (marked with the blue diamond), and not 60 s as one would expect.Clearly, the ramp rate is not linear.Furthermore, let us assume that this turbine-generator participates in FCR provision: in the first case one could estimate the ramp rate is 5 kW/30 s = 10 kW/min, while it the second case it is 10 kW/35 s = 17 kW/min.Thus, by using the first ramp-rate one could deduce that the generator cannot provide more than 10 kW/min × 0.5 min = 5 kW reserve while in reality it could probably provide more.Similarly, by using the second ramp-rate, one could deduce that the generator cannot provide more than 17 kW/min × 0.5 min = 8.5 kW reserve, while it clearly can because at T = 30 s the output power is ≈ 11 kW which is greater than 8.5 kW.By using the second ramp-rate, one could also deduce that by offering 5 kW of reserve, the limit will be reached in proportionally smaller time 5 kW ÷ 17 kW/min = 17.6 s, while it clearly takes around 30 s no matter the size of the power set-point change.Therefore, our main argument is that when the time constant of a generator is of the same order of magnitude as the FCR activation time, constant ramp rate assumption may lead to erroneous conclusions that a generator can provide FCR while in reality it cannot, or vice versa.That is why we take the approach of using the settling time instead.
The dynamic behaviour of units is nonlinear due to delays and time lags of its components [20].Without loss of generality, this behaviour in most cases can be described by the first-order or Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the second-order transfer functions (i.e.exponential behaviour or damped sinusoidal behaviour, respectively).This means that the rate-of-change-of-power is not constant nor linear.The rate limiters only constrain the maximum opening and closing speed of the gates or valves of conventional generators, or constrain the actuation speed of power electronic converters.However, these slew rates are usually much greater than the time constants of the whole system and will not be a constraining factor during normal operation.Thus, the unit behaviour will be determined by the slowest time constants.For example, opening the valves of a steam turbine very quickly is useless if the time it takes to translate the valve opening into power at the turbine shaft is very slow (due to slowness of the thermodynamic processes).The only way to achieve a linear ramp behaviour is to set these rate limits much slower than the inverse of the unit's largest time constant, but this is usually not done for reasons related to performance and safety, especially in the context of providing FCR.Hence, most grid codes do not require the change of output power of FCR-providing units to be linear during the activation period, but only that the full FCR volume is activated within a prescribed time [2], [16], [17], [21], [22], [23].For a more detailed discussion on this topic with examples, we refer the reader to [24].
To implement this nonlinear dynamic behaviour of FCRproviding units into an aggregator's MILP model for bidding in the FCR and the DA markets and to ensure grid compliance during the FCR activation, we propose incorporating the settling time constraints of FCR-providing units.Settling time is the time required for a unit to reach a new set-point within a certain accuracy (99.5% in this article).Settling time constraints can be derived analytically if a dynamic model is known or they can be experimentally obtained on-site from simple step response tests.

C. Dynamic Model Formulation
We derive settling times of hydro, thermal and converterinterfaced storage units from their commonly used simplified dynamic models.These settling times are then used as constraints in the MILP model in Section II-D.Although the used dynamic models are relatively simple, they are adequate for this type of study, as proven in [18], [25], as they include all relevant nonlinearities per grid-code requirements [2]: the insensitivity, the deadzone, the rate limits and the saturation.The presented procedure is valid even if more complex dynamic models are used.
For better clarity, indices u, d, t, w that are used in the optimisation model (13a)- (30) are omitted in the dynamic model ( 1)- (12).
1) Hydro Turbine-Generator Portfolio Model: The aggregator's HTG portfolio consists of N HTG hydraulic turbines with a digital PID governor modelled by the equivalent transfer functions shown in Fig. 3.The model also considers the gate opening rate limits and the saturation, as well as the control system insensitivity and the deadzone.p DA u is a power set-point of hydro unit u ∈ {1. ..N HTG }, which is modified in time according to a determined day-ahead plan.The model from Fig. 3 can be used to model various kinds of hydraulic turbines of different sizes.More info on hydraulic turbine-governor modelling can be found in [25], [26], [27], [28], while in this work the details are omitted for brevity.
It suffices to say that T w is the water starting time determined by the hydraulic turbine design as well as the operating point.Time constants of the equivalent governor transfer function T A , T B and T C depend on the chosen control system parameters.For a digital PID hydraulic governor, these time constants are defined as follows: where K p , K i and K d are proportional, integral and derivative gains of the PID controller; T r is the washout filter time constant (equivalent to the dashpot/reset time constant of classic mechanical-hydraulic governors); k is the droop gain of a unit and T g is the governor time constant (main servo).The time constants of the turbine governor will be greater than time constant T w of the turbine for a standard range of parameters, thus the response of the hydro turbine-governor will be determined mostly by the governor [26].If the turbine dynamics are neglected, the behaviour of a single hydro generator can be approximated by a first-order transfer function ( 4), assuming T B T C [25].
For a step change of the input Δu(s) (Fig. 3), the settling time of a new set-point is equal to Taking into account the induced errors due to simplification of the model, a safety margin of 10 time constants instead of 5.3 can be used to guarantee a hydro unit will reach a new set-point within ±0.5% accuracy in T HTG : Therefore, to decide whether a specific hydro unit can participate in the FCR provision, it is sufficient to check whether T HTG = 10 T B T A ≤ T FCR in the linear programming framework, where T FCR is the FCR activation time, usually 30 seconds.Since T HTG depends on the droop gain per (1) and (2), and droop gain is also a decision variable in the linear programming framework, (5) can be rewritten as a function of k using (1) and (2): Based on (6), a constraint for HTG droop gain k (in MW/Hz) can be derived: where P max /f n is used to convert from per-unit to MW/Hz.Equation ( 7) corresponds to constraint (21) in the MILP model.

2) Thermal Turbine-Generator Portfolio Model:
The aggregator's STG portfolio consists of N STG steam reheat turbines with a hydraulic governor modelled by equivalent transfer functions shown in Fig. 4. The model also considers the valve opening rate limits and the saturation, as well as the control system's insensitivity and the deadzone.p DA u is a power set-point of thermal unit u ∈ {1. ..N STG }, which alters in time according to a determined day-ahead plan.This model is suitable for describing e.g.small, distributed biomass generators.More info on the steam turbine-governor modelling can be found in [25].The dynamics of the governor in thermal generators are much faster than the dynamics of the turbine with the largest time constants due to the slowness of the thermodynamic phenomena.Thus, if the governor dynamics are neglected, the dynamic response of a thermal generator can be approximated by: where Δu is the power set-point change due to frequency droop control (Fig. 4), T R is the reheat time constant and F H is the portion of the turbine power developed at the high-pressure stage.Settling time of the new set-point can be derived exactly and is equal to (9).For a thermal generator to be eligible to participate in the FCR, T STG has to be less or equal to T FCR (equations ( 22) and ( 30)).Unlike hydro units (see eqs. ( 6)), in (9) the settling time of a thermal unit is independent of the droop gain k: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.3) Energy Storage Portfolio Model: Batteries and supercapacitors are static energy storage systems interfaced to the grid via voltage source converters.There is no mechanical motion in those devices, thus the response is very fast and determined by the time constant of the inverter control system.Behaviour of these devices can be modelled as a first-order transfer function [29], [30], [31].The same model is used for modelling both the battery and the supercapacitor portfolios and is provided in Fig. 5.It considers the programmed rate limits and the saturation, as well as the control system's insensitivity and the deadzone.p DA u is the power set-point of storage unit u ∈ {1. ..N BSS/SCS }, which is modified in time according to a determined day-ahead plan.
The charge/discharge logic prevents charging/discharging when the SOE is at the maximum/minimum.The round-trip efficiency is used and is assigned to the charging mode only [32].It is modelled with efficiency function η(p u ), where p u is the output power of storage unit u: The dynamic response of a BSS/SCS can be approximated by: where Δu is the power set-point change due to frequency droop control (Fig. 5) T c is the converter time constant.Consequently, the settling time of a new set-point is equal to (12), and T ESS ≤ T FCR is a prerequisite for the FCR provision.The settling time of a static storage unit is also independent of the droop gain k.Since the time constants of these devices are in the order of several hundred milliseconds, this will not be a constraining factor.Additionally, if a linear response is desired, slower rate limits can be used while still providing fast service.( 12) corresponds to eqs. ( 22) and ( 30) in the MILP model.

D. Mixed-Integer Linear Programming Formulation
The objective of the MILP model is to maximize the aggregator's profit from trading in the FCR and the DA markets.The first term in (13) represents the income from the FCR market, the second one represents the income from the DA market, the third term (13c) represents compensation for the activated FCR energy, remunerated at the DA price, while the last term (13d) is the cost for recovering the BSS SOE at the DA price due to its round-trip inefficiency.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

−
h∈H t|t∈h u ∈U w∈W subject to: Set U = {HTG, STG, BSS, SCS} is a set of DER technologies at the aggregator's disposal.Sets H and T represent timesteps, T are hours in a day, while H are six 4-hour intervals for FCR. 1 Therefore, t ∈ h denotes hours belonging to each interval h.D is a set of representative days, representing characteristic features of different FCR and DA market price correlations, and W is a set of scenarios of the FCR activation on the day of delivery.
The aggregator's decisions in the FCR market are modelled in constraints ( 14)-( 16).( 14) calculates the reserved power capacity of the aggregator in every time-step by multiplying the droop gain of each unit of each technology by the frequency deviation.The frequency deviation is defined as Δf = Δf max − dz, where Δf max is the maximum frequency deviation (at which FCR is fully activated) and dz is the combined effect of deadzone and insensitivity.Inequality (15) imposes that the reserved capacity does not exceed the 30-second ramping potential of the units that provide the FCR.Binary variable y ut is equal to 1 if a unit is used for FCR and to 0 if not.Constraint ( 16) locks the droop gain of each unit to a single value during each 4-hour period.
All units, regardless of the technology, are subject to constraints ( 17)-( 21).( 17) and ( 18) constrain the power output of each unit by its minimum and maximum capacities if the unit is on, which is determined by binary variable x ut .For the storage systems, this variable can be omitted because the value of P min u is negative, i.e.P min u = −P max u .Furthermore, these two constraints do not contain the first term (p FCR udtw ) for SCS since they do not participate in the DA market due to their low available energy.Inequalities (19) and (20) constrain each unit's droop gain between its maximum and minimum values.Constraint (21) defines the maximum droop gain setting K T based on a HTG unit's settling time dynamics from eqs. (6).Constraint (22) prevents the STG, BSS and SCS units, having settling times longer than 30 seconds, from offering the FCR through binary variable y ut , with Y u defined as follows: Energy storage technologies have additional constraints ( 23)-( 29).( 23) defines variable p DA ut as a difference between the discharging and charging powers.For SCS, p DA ut , ch DA ut and dis DA ut are zero, again due to their low available energy.Constraints ( 24) and ( 25) ensure that storage systems do not charge and discharge simultaneously.( 26) calculates the SOE of storage technologies considering the SOE in the previous time-step, the DA charging and discharging powers, the FCR activation energy Δe FCR udtw and the SOE recovery for FCR market defined as the opposite value of SOE change from the previous time step Δe FCR ud,t−1,w .Constraints ( 27) and ( 28) limit the SOE of storage systems considering the maximum FCR activation.Finally, the change in SOE resulting from the actual FCR activation is calculated in (29) by multiplying the droop gain k udtw , maximum frequency deviation Δf , time step duration Δt and parameter âadtw which represents percentage of FCR activation for each unit.

III. STUDY CASES
To demonstrate performance of the model described in Section II, we study a case of an aggregator operating a heterogeneous set of DERs.We consider forty units in total with randomly-generated parameters divided into four technologies: 10 × HTG, 10 × STG, 10 × BSS, and 10 × SCS.The installed capacity of all units ranges between 0.1 MW and 10 MW, totalling 44.75 MW of HTG, 52.9 MW of STG, 47.1 MW of BSS and 57.3 MW of SCS.The energy-to-power ratio for BSS is 1, while SCS have discharging times between 5 s and 30 s. Bounds for the units' droop settings are 2-6% for hydro units, 2-7% for thermal units, and 0.4-12% for batteries and supercapacitors.The 40 units operate in a power system of a size comparable to the interconnection of continental Europe (≈ 1000 GW of installed generation capacity and consumption in the range 300-600 GW [34].Nominal frequency for continental Europe, which is considered in this article, is 50 Hz, with the maximum steady-state frequency deviation 200 mHz, the combined effect Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 6.Modified IEEE 14-bus system. of deadzone and insensitivity of 8 mHz (which is under the permissible limit defined by [16]).
The largest European coordinated FCR market has a product delivery period of 4 hours, traded on the day prior to the delivery, the products are symmetrical, and the minimum bid and product resolution are both 1 MW [33].The energy-only DA market usually has a 1-hour resolution and is traded on the day prior to the delivery, same as the FCR market.This market setup is implemented in this article.The resulting optimal portfolios are tested for the ability to provide FCR according to the technical requirements set by Tennet [2].
To do this, we use a modified IEEE 14-bus system to emulate the continental Europe synchronous area in terms of its size.All the generators are equipped with turbine governing systems and automatic voltage regulators.We expanded the IEEE 14-bus system (Fig. 6) with a 35 kV level distribution feeder radially expanded from Bus 12 where our portfolio of 40 units is distributed (Fig. 7).Thermal and hydro units in our portfolio are equipped with corresponding turbine-governing systems (TGOV1 and HYGOV, which correspond to simplified models in Figs. 3 and 4) and excitation systems (IEEET1S).The battery and supercapacitor systems are modelled as gridfollowing converter-interfaced units that include outer active (Fig. 5) and reactive power control loops, inner current control loops as well as phase-locked loop dynamics for grid synchronisation [29].The seventh order dynamic model is used for all the synchronous generators in the system.The dynamic model of the transmission and distribution systems is developed in DIgSILENT PowerFactory and simulated using the fundamental frequency (RMS) simulations with a 0.1 ms time step.Main parameters of the 40 aggregator's units are provided in [35].
The developed models are used to examine the motivation behind the aggregator's behaviour with regard to different parameters of the aggregated DERs and to analyse the performance of such entity in a real-time setting.We assume the aggregator optimises the droop gain setting for the units on a day-ahead basis and signals the changes in settings on the 4-hourly or the 1-hourly basis.

A. Description of Study Cases
The initial step in the process illustrated in Fig. 1 was the standard scheduling model considering the reserve capacity allocation, but without activation (see Table I).The storage units scheduled using this model were not able to deliver the FCR in full because they did not have enough stored energy [2].For this reason, constraints ( 27) and ( 28) were expanded to their current form, ensuring that at the start of each interval there is enough energy stored for delivering the maximum FCR activation in both directions.The model also lacked ramping and settling time constraints for generating units offering FCR, so any unit, even those not fast enough to deliver FCR in time, were scheduled to provide FCR.This was corrected by adding (15), which is a standard ramping constraint used in unit commitment problems.However, this constraint was insufficient to ensure compliance with the grid code, so constraints ( 21) and (22) were added.After reaching the grid compliance, we consider the following five model formulations based on the possible real-world FCR market setups: 1) F -The base case where the FCR market participation, i.e. selling FCR capacity in the market, is modelled without considering the activation and its recovery, i.e.SOE variations during the real-time operation.This model is used in most studies in the literature that consider FCR market participation.While such model might be correct when considering only generating resources providing FCR, it might not be accurate for models including energy storage, where the FCR activation is omitted due to the assumption that the upward and downward FCR activations cancel out over the time.However, even if that was the case, energy storage with non-ideal cycle efficiency would deplete with time.2) FA -The FCR activation is considered in the SOE calculation constraint, but not remunerated in the objective function.This case is used in several recent studies on the FCR market-participating energy storage, including [9].The SOE recovery is not considered.3) FAR -The FCR activation is treated the same as in the previous case, with the addition of the required SOE recovery.This case is arguably the most realistic setup under the current market rules.4) FAD -The activation of FCR is compensated at the DA market price.The activation is included in the model in the same way as in the FA case, i.e. with remuneration included in the objective function.The SOE recovery is ignored.5) FADR -The FCR activation and its compensation is treated the same as in the previous case, while the SOE recovery is treated as in the FAR case.Cases 3) -5) were not used in the known literature.Cases 4) and 5) are the furthest from the current state-of-the-art, as the FCR activation is currently not compensated in European markets.We use the cases to point out the differences between the FCR modelling approaches.Table I provides a quick reference of the constraints used in each of the five cases.

B. Data
We used the French FCR market data for year 2020.The reserved capacities and activated energies were used to generate percentages of hourly FCR activation.For the DA market prices we used the Croatian Power Exchange (CROPEX) data for year 2020.We prepared data for the study cases using the k-means method described in [36], clustering the data into four representative days considering the DA prices, the FCR prices and the capacities.The optimal number of representative days, four, was determined using the elbow method.
We further subdivided each representative day into five scenarios considering the upward and downward FCR activations.Activation of FCR was assigned to each unit to capture the effects of energy storage efficiency.The FCR activation parameter was calculated by dividing the known parameters: activated FCR energy and reserved FCR capacity over an hour.To create sub-scenarios, we followed the same procedure with FCR activation as with the DA and FCR capacity market prices, using the k-means clustering algorithm with the elbow method for optimal number of clusters.Fig. 8 represents the four scenarios and corresponding sub-scenarios.

A. Dynamic Model Validation
First, we show that the simplified constraints introduced in Section II-C sufficiently accurately describe nonlinear unit dynamics in the time horizon of interest.Therefore, the simplified models in Simulink ((4), ( 8) and (11)) are compared to both nonlinear models in Simulink (Figs. 3-5) that consider only active power dynamics as well as to full nonlinear model in DIgSILENT PowerFactory (Figs. 6 and 7) that includes grid topology and high-order device models.
Droop gains were randomly assigned to all 40 units in the aggregator's portfolio.A −0.2 Hz step change was applied to units' active power control systems in order to emulate one of the main prequalification tests for FCR provision [2].Results are shown in Fig. 9.It can be seen that the introduced simplifications (dashed orange lines) cause inaccuracy only during the initial transient which is permissible since our methodology only looks at the settling time, i.e. whether the required power will be delivered in under 30 s activation time.Simplified models behave identically to nonlinear models after 10 s for steam turbines, after 20 s for hydro turbines and after 1 s for static energy storage units which is below the required activation time.Therefore, the derived simplified models are adequate.

B. Prequalification Procedure
In this section we show that the devices selected by the MILP framework ( 13)- (29) for FCR provision indeed satisfy the prequalification requirements imposed by Tennet [2].There Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.are several tests that need to be performed sequentially, but due to limited space only a response to a step change in frequency deviation from 0 to −200 mHz is shown here, as illustrated in Fig. 10.Each test was performed for all cases listed in Table I and each model was compliant to the requirements.For instance, if the settling time constraints were calculated incorrectly, this would show in the results of the dynamic simulation.Since the settling times of the units were approximated by neglecting the nonlinearities and all but the largest time constants, it is possible that the error induced by the approximation makes the FCR activation non-compliant (this depends on the used model and the order of approximation).For the parameters used in this article there were no non-compliant cases.If there were, that would be an indicator to re-tune the controller parameters (if possible) or increase the safety margin on the settling time constraints.
Let us consider an aggregator who wants to provide 10 MW of FCR.If the settling time constraints ( 21)- (22) are not included in the optimization framework, then a portion of these 10 MW FCR may be assigned to units that are too slow to provide it regardless of their rate limits.This would result in delivering less FCR than contracted, which would incur penalties for the aggregator.At the FCR activation stage, the reserved capacity would not be completely activated within T FCR = 30 s, as illustrated by the solid blue line in Fig. 10.With the settling time constraints included, the optimization framework dispatches the FCR capacity only between units capable of delivering the full capacity within T FCR .Therefore, all cases listed in Table I are constrained by ( 21)- (22).The same process was repeated for constraints ( 23)- (29).
Note that the MILP formulation with settling time constraints (21)-( 22) is conservative.It discards any unit with settling time greater than the activation time T FCR from providing a share of full pre-qualified FCR volume.However, it is possible that the aggregated dynamic response of the FCR-providing units could satisfy the required performance within T FCR even though some units have settling times longer than T FCR .However, due to the limited space, a more detailed investigation is omitted in this work.

C. Scheduling
In this section we present the results of the MILP model.Fig. 11 shows the FCR and DA schedules of the full portfolio, as well as the BSS SOE for Case F, while Fig. 12 provides the same information for Case FA.The first conclusion from these figures is that the SCS are not scheduled to offer FCR at any point.In fact, this is true for all five cases listed in Table I.This outcome is explained by the FCR market rules, which require that an energy storage unit has sufficient energy stored to provide reserve for frequency deviation of ±50 mHz for 15 minutes, followed by 15 minutes at full deviation capacity [2], while SCS at full power can provide reserve for up to 30 s.Following the FCR provision rules, this amounts to a minimum of 31.25% of p FCR h •ΔT .The SCS energy-to-power ratio is by their physical properties very low, which makes the FCR market participation infeasible for this technology.
As all units have finite capacities and the units' production costs are ignored, the division of the capacity between the two markets can be seen as a zero-sum game.The generators are scheduled to participate in the FCR market with their maximum possible capacity constrained by their droop gain limits, while the rest of the capacity can be scheduled to the DA market.While schedules of the HTG units do not vary throughout the day, the STG and BSS schedules do.Even in Case F, where FCR activation is not included (Fig. 11), the BSS take an opportunity to benefit from price arbitrage in the DA market, while offering most of their capacity for FCR provision.This is evident from Fig. 11(c), which shows the BSS are charging in hours 8-9, while in the next 4-hour period the FCR bid does not increase.We can therefore conclude that the arbitrage is the sole purpose of the BSS charging in this period.The variation in the BSS schedule is more pronounced for Case FA, shown in Fig. 12, with activation of the FCR included in the model.Fig. 12(c)) illustrates how, as the stored energy is impacted by the FCR activation, the BSS show more variation in the FCR bids and the BSS SOE limitations become more prominent.This causes the BSS to take a more active role in the DA market which, in addition to arbitrage, they use for SOE recovery.

D. Comparison of Study Cases
Economic comparison of the cases is presented in Table II.The case without dynamic constraints and without the FCR activation (bottom row) is taken as a reference, as it represents the state-of-the-art approach.Below each objective value is the percentage change in relation to the reference case.Our proposed modelling approach (row three) achieves lower profit for each of the cases.The reason for this is the lack of dynamic constraints in the standard models, which allows the non-grid-compliant generators to provide FCR, increasing the aggregator's profit.Comparison to the DA-only and the FCR-only cases makes it evident that most of the profits are generated in the DA market.This is not surprising, since FCR prices are around 10% of the DA prices.Another positive outcome of our approach, compared to the standard approach is that solving the model with dynamic constraints takes 4%-24% less time, depending on the model.
For the synchronous generators, the DA market participation is a prerequisite for the FCR market because a unit must be online and generating electricity to offer reserve.For storage units this is not the case, as can be seen in Fig. 13, which provides further comparison of the cases.
Each dot in the plot represents one instance of the model results.The total number of dots corresponds to 4 representative days, 24 hours, 5 scenarios and 10 units, totalling 4800 points.The color of each dot represents the number of times the model resulted in that combination of the DA and FCR market bids.The generators' FCR capacity bids are constrained by their ramping capabilities based on the rate limits and the settling time constraints on the one hand and by their droop gain limits on the other.Out of the ten STG in the aggregator's portfolio, by assigning random parameters, three ended up having the settling time that allows for the FCR market participation.The top row of graphs in Fig. 13 represents the HTG.Their bids are also constrained by their droop gain limits, which translate to the FCR capacity.This reflects on the units' schedules so the HTG and STG offer at most 20% of their power in the FCR market.This percentage comes from the minimum droop gain coefficient, which is 2% for both technologies, so their maximum droop gain is K max = 0.192 p.u./Hz, where the maximum frequency deviation Δf is 192 mHz after the 8 mHz combined effect of deadzone and insensitivity is removed from the 200 mHz required by the Regulation EC 2017/1485 [16].Their remaining capacity is offered in the DA market because the generating units need to be online (symmetrical provision of FCR is assumed) to provide FCR.The most common decision for HTG is to offer 100% in the DA market.For generators, the amount scheduled to the FCR market is largely dependent on the units' production costs, which is why the cheaper units (HTG, top row of Fig. 13) schedule all available power to the DA market, even though their settling times allow the FCR market participation.On the other hand, STG have higher production costs and therefore the aggregator finds it profitable to divide their capacity between the two markets.
The bottom row of graphs in Fig. 13 represents the BSS, whose droop gain limits are set lower so that they do not affect their FCR market bids, allowing them to offer their full capacity for FCR provision.Moreover, their bid FCR capacity increases as the DA volume decreases.Thus, the FCR reaches the maximum when the DA is equal to zero.At that operating point, they are constrained only by the stored energy.Fig. 13 shows how the five proposed cases listed in Table I can be divided into two groups with regard to storage units.The first group are Cases F, FAR, and FADR, while the second one are Cases FA and FAD.In the first group, the FCR market bids are concentrated around the 1 p.u. point, while in the second group the FCR bids are more evenly spread between 0 and 1 p.u.The first group comprises cases where the FCR activation is not considered or the activated energy is recovered within the next hour.The second group contains the cases where the FCR activation is considered, but SOE recovery is not.The models in the first group provide similar results with respect to the distribution of capacity between the FCR and the DA markets, with those including the FCR activation exhibiting more variability in the FCR and the DA bids, indicating that they can be used interchangeably.The F model, without the FCR activation, would be useful in applications where hour-to-hour SOE fluctuation is not relevant, while the FAR should be useful in operational models where the battery's SOE has more impact on the schedule.More variability in the BSS schedules can be observed in the second group, i.e. the cases with activation but without energy recovery.These models demonstrate more risk-averse behaviour of the BSS in terms of the FCR market bids, caused by the FCR activation uncertainty and inability to recover SOE through trading on the day of delivery.In terms of the DA market bids, these models explicitly schedule more power from the DA market for the BSS charging, which is in the first group of cases included in the model on an assumption that trading is performed on the day of delivery, i.e. in the intraday market.

V. CONCLUSION
The aggregator in this article centralises the scheduling of the DER units under its control and bids on their behalf in the FCR and DA markets.The risks of trading in these markets are captured by stochastic programming based on three sources of uncertainty: the FCR capacity prices, the DA market prices and the FCR activation.
The goal of this article was to harmonise the decisions made at the operational level with the dynamic limitations of the aggregated units.Therefore, we presented a methodology for using a detailed dynamic system framework to incorporate the dynamic constraints within the MILP model.This ensured the ability of the aggregator's portfolio to deliver all contracted services to the system.Using the dynamic system model we demonstrated that the aggregator's portfolio complies with the prequalification requirements set by the system operators in all considered cases.
The results illustrate that including the FCR activation in operational models shifts more capacity of storage technologies from reserve to energy markets, which can be seen as a risk mitigation measure.However, not all modelling approaches gave the same results when the FCR activation was included.The models that include the DA and FCR market participation as well as reserve activation, but ignored the need for SOE recovery, resulted in larger volumes traded in the DA market.This demonstrates a more risk-averse behaviour in comparison with the cases that did not include the SOE recovery.In these cases, the aggregator charges or discharges the BSS in advance to be able to provide the correct amount of FCR service, i.e. the SOE recovery is performed in advance.
generator.Manuscript received 18 September 2022; revised 15 March 2023; accepted 2 July 2023.Date of publication 5 July 2023; date of current version 21 February 2024.This work was supported in part by the European Union through the European Regional Development Fund for the Competitiveness and Cohesion Operational Programme 2014-2020 of the Republic of Croatia under Project Number KK.01.1.1.07:Universal Communication and Control System for Industrial Facilities and the Project evelopment of a System for Control of Electrical Energy Consumption in Households SUPEER, under Contract Number KK.01.2.1.02.0063.The work of Ivan Pavić was supported in part by the Luxembourg National Research Fund (FNR) and PayPal, PEARL grant Reference 13342933/Gilbert Fridgen and in part by the Luxembourg National Research Fund (FNR) -FlexBeAM Project, Reference 17742284.Paper no. TPWRS-01382-2022.(Corresponding author: M. Miletić.)

Fig. 1 .
Fig. 1.Proposed methodology for ensuring grid-code compliance of the MILP results.

Fig. 2 .
Fig. 2. Step change of active power for a laboratory-scale hydro turbinegenerator.

Fig. 8 .
Fig. 8.The four DA and FCR price scenarios with their corresponding FCR activation sub-scenarios.

Fig. 10 .
Fig. 10.One of the prequalification tests with and without the settling time constraint (21) and (22).

Fig. 11 .
Fig. 11.Expected DA and FCR bids, as well as BSS SOE, for Case F.

Fig. 12 .
Fig. 12. Expected DA and FCR bids, as well as BSS SOE, for Case FA.

Fig. 13 .
Fig. 13.DA vs. FCR market bids for HTG and BSS (BSS DA market bids range from -1 for maximum charging to 1 for maximum discharging power).
Optimal Primary Frequency Reserve Provision by an Aggregator Considering Nonlinear Unit Dynamics M. Miletić , Student Member, IEEE, M. Krpan , Student Member, IEEE, I. Pavić , Student Member, IEEE, H. Pandžić , Senior Member, IEEE, and I. Kuzle , Senior Member, IEEE

TABLE I CONSTRAINTS
USED IN THE MODEL BY CASES

TABLE II COMPARISON
OF PROFITS BETWEEN CASES IN EUR