Monetarization of the Feasible Operation Region of Active Distribution Grids Based on a Cost-Optimal Flexibility Disaggregation

Hierarchical grid control strategies are an appropriate design concept for the coordination of future transmission system operator (TSO) and distribution system operator (DSO) interactions. Hierarchical approaches are based on the aggregation of decentralized ancillary service potentials, represented by converter-coupled, communicable active and reactive power flexibility providing units (FPU, e.g. wind turbines) at vertical TSO/DSO system interfaces. The resulting PQ-polygon made available by the DSO for a potential request of ancillary service flexibilities by the TSO is called feasible operation region (FOR). A monetarization of the FOR is necessary for the implementation as operational degree of freedom within TSO grid control. In the context of a local DSO and global TSO market, this article presents an approach for the monetarization of the FOR by a cost structure using metadata from population based aggregation methods. At the local DSO market free bids for the active and reactive power flexibilities by the FPUs stakeholders are assumed. Within the aggregation method multiple FPU flexibility polygons at a single bus are aggregated for a reduction of the search space dimensions. Thereby, the main contribution of the proposed method is the cost-optimal disaggregation of a flexibility demand to the single FPUs within the aggregated FPU by a mixed integer linear program. The approach can be simply adapted for the coordination of DSO/DSO-interactions regarding hierarchical multi-level grid control strategies.


I. INTRODUCTION
The transition of the electric energy system leads to a massive integration of decentral energy resources (DER), especially to the distribution system level [1]. At the same time the contribution of conventional thermal power plants to the energy mix is decreasing in Germany, due to the coal and nuclear power phase-out and the priority feed-in of volatile renewables [1]. The ancillary service potentials of the transmission system operators (TSO), guaranteeing a safe and reliable energy supply, are reduced. Additionally, there is a need for more control measures (e.g. active and reactive power redispatch) to avoid grid congestions and voltage band violations resulting from a local power imbalance [1]. A variety of system elements at the distribution system level are converter coupled and have an information and The associate editor coordinating the review of this manuscript and approving it for publication was Ravindra Singh. communication interface [2]- [7]. Thereby, they are flexibly controllable within the system operation of the corresponding distribution system operator (DSO). The possible adaptation of the active and reactive power supply of these flexibility providing units (FPU) can be described as polygon at the PQ-plane (see Fig. 1) [2], [3], [8]. The distribution system level transforms and becomes increasingly active. The flexibilities within active distribution networks (ADN) can be used for a change of the vertical active and reactive interconnection power flows (IPF) at the TSO/DSO system interface according to the demand of the TSO and by this as additional ancillary service potential [9]- [14]. To coordinate the technical and organizational interactions between the system levels an amendment of the regulatory framework and the bilateral agreement between TSO and DSO is necessary [5], [7], [12], [15]. In literature especially hierarchical grid control strategies are discussed for the coordination of future TSO/DSO-interactions [16]- [18]. In general, these approaches are based on a dayahead or intraday, proactive determination of a feasible operation region (FOR) (see Fig. 1). Therefore, the DSO aggregates the flexibility potentials of the flexibility providing units connected to the distribution system level at the TSO/DSO system interface (see I. in Fig. 1) [2], [9], [14], [19]. The FOR describes the possible adaptation of the active and reactive IPFs (P vert , Q vert ) by the DSO within a PQ-plane. The next step is the specification of the flexibility demand (P vert,sp , Q vert,sp ) of the TSO from the distribution system level within the operational management of the TSO (see II. in Fig. 1). Finally, the TSO flexibility demand is distributed to the FPUs within the operational management of the DSO (see III. in Fig. 1). Thereby, the TSO/DSO-coordination process is based on bottom-up aggregation, top-down specification and local distribution and can be also applied for the coordination of DSO/DSO-interactions [2], [3], [16]. Thereby, a cascading process for the coordination of system operator interactions within the overall system results in the context of multi-(voltage-)level grid control strategies [16].

II. STATE OF THE ART OF FOR MONETARIZATION
The determination of the FOR at a single vertical TSO/DSO system interface is based on an aggregation of the individual flexibility areas of the FPUs connected to the DSO system. Rudimental approaches estimate the FOR (see Fig. 2a) or use the Minkowski sum (see [20], [21]) for the aggregation of polygonal flexibility areas of the FPUs, neglecting security constraints (e.g. voltage limits, maximum thermal currents) [22].
In literature, two main categories of advanced aggregation methods also considering security constraints, the active and reactive power demand of the DSO system as well as variations of the voltage at the interconnection bus i as three-dimensional PQV-FOR (see Fig. 2b) are currently investigated. These are stochastic approaches using random search (RS) and optimization based approaches solving an adapted optimal power flow (OPF) problem (cf. [2], [23]).
RS aggregation methods (cf. [2], [14], [19], [24]- [26]) compute a variety of power flow calculations (e.g. Newton Raphson) for different constellations of operating points of  the FPUs as Monte-Carlo scenarios. Each solution of a Monte-Carlo scenario represents a point at the PQ-plane which results in a point cloud. Disadvantages of RS methods are the long computation time and a challenging determination of the specific FOR edge especially in cases of nonconvexities [14]. Straight forward RS approaches neglecting the covariance between the FPUs lead to a convolution of the results based on the central limit theorem of the probability theory [24], [27]. Novel RS approaches are using probability density functions for an appropriate sampling of the FOR [2], [28].
Advantages of RS approaches are the generation of metadata (e.g. flexibility provision per bus, bus voltages) for each point within the FOR (see Fig. 3a) and the simple consideration of non-convex FPU flexibility polygons [2], [23]. Another advantage is the good performance of the RS approaches for larger systems. The individual computation time for one power flow calculation increases proportional to the number of buses [2], [19]. Thereby, the computation time scales predominantly by the number of power flow calculations.
OPF-based aggregation methods (cf. [9]- [11], [28]- [32]) sample the FOR edge in incremental steps. Therefore, at each sampling step an OPF with a specific constellation of the active and reactive IPFs is solved. The significance of the advantages of the OPF-based approaches depends on the specific optimization method for the solution of the OPF problem. Linear or quadratically constrained linear programming approaches (LP) determine the OPF edge in short computation times [8], [9], [23], [30], [33], [34]. The non-linear system behavior is only approximated which can lead to local VOLUME 10, 2022 convergence and by this to an over-or underestimation of the FOR compared to the FOR resulting from other optimization methods [2], [8]. This over-and underestimation of the FOR needs to be minimized to provide the TSO as much guaranteed flexibilities as possible in case of critical system states. The performance of LP scales proportionally for larger systems but the quality of the results decrease due to the increasing number of non-linearities [30], [35], [36]. Nonlinear programming (NLP) approaches lead to an appropriate sampling of the FOR in short computation times also identifying non-convexities of the FOR [9]- [11], [23], [25]. The performance and the quality of the results of the NLP depends on the system size and the specific solver (e.g. interior point optimizer). With an increasing number of operational degrees of freedom and constraints also the possibility of local convergence increases [11].
The information of the cost structure of the FOR can be used by the TSO for an economical specification of the required flexibilities at the vertical system interconnection to the DSO system within a local DSO and global TSO flexibility market [9], [14], [37]- [39]. RS as well as OPF-based aggregation methods are compared with each other in literature (see [9], [14]) for the determination of the cost structure of the FOR. The description of the cost structure depends on the individual aggregation method. For RS approaches each point of the resulting point cloud at the PQ-plane can be simply monetarized based on the metadata received by the multiple power flow calculations (cf. [9]). The interpretation of the cost structure and the identification of zones with a uniform price is challenging (see Fig. 3a). The reason for this are multiple overlapping points resulting from various possible constellations of FPU flexibility provisions for a specific IPF. The resulting cost structure represents an estimation for the costs that the TSO can expect for a specific active and reactive power demand from the DSO system.
Solely the prices of the FOR-edge are available after the sampling process for OPF-based approaches in contrast to RS approaches (cf. [14]). Multiple sampling processes that consider the compliance of specified prices as additional constraints are necessary to determine the cost structure within the FOR [14]. The cost structure is described by price zones z whose individual contour lines represent a specific price c max,z (see Fig. 3b) [14]. An advantage regarding RS approaches is the guarantee of minimum costs for the contour lines. Disadvantages are that the cost structure between two contour lines is not available. The computation time of OPF-based approaches is increased for an appropriate sampling of the cost structure.
Approaches using a metaheuristic (e.g. particle swarm optimization, genetic algorithm) combine the advantages of RS in generating metadata and the appropriate sampling of the FOR-edge of OPF-based methods [9], [23], [28]. In general, solving an OPF problem by metaheuristics is based on an evaluation of the objective function by a variety of power flow calculation, analogously to RS approaches. An iterative, algorithm-specific adaptation of the population is used for convergence in a solution of the OPF problem.
As typical for stochastic approaches, metaheuristics cannot guarantee global convergence or determine the remaining gap to the global optimum [40]. Nevertheless, the investigations in [23] show a good performance of the particle swarm optimization (PSO) in determining the edge of the FOR and a high quality of the results (e.g. size of the FOR, identification of non-convexities) compared to NLP and quadratic constrained LP.
This article continues research concerning an application of stochastic and metaheuristic aggregation methods in sampling the FOR at the TSO/DSO system interface using the PSO. The focus is on the determination of the FOR cost structure within a local DSO and global TSO flexibility market using metadata generated by the PSO [18], [37], [41]. The main contribution is the consideration of multiple FPUs that are aggregated by the Minkowski sum at a single bus with individual active and reactive power cost bids. First, the flexibility potentials of FPUs and the general process of the FOR determination by the PSO are described in section III. The method is introduced for non-convex FPU PQ-polygons and for the example of a two-dimensional FOR. The process can be simply adapted for the determination of a threedimensional PQV-FOR due to slack voltage variations. The extensions regarding the FOR monetarization are described in section IV as a mixed integer LP (MILP), to specify the costoptimal distribution of a flexibility demand to the individual FPUs per bus. The MILP can be reduced to LP in case of only convex FPU flexibility polygons. The general principles of a local DSO and global TSO flexibility market based on commodity and service prices for FPU flexibility provision are introduced for the monetarization of the FOR. The presented process can be simply adapted for other population based aggregation methods like RS or further metaheuristics. Furthermore, the process is applicable for the coordination of DSO/DSO-interactions. Therefore, the case studies in section V are based on MathWorks MATLAB simulations using an adopted Cigré MV test system [42]. The dataset of the grid model including the FPU PQ-polygons is accessible at [43] for reproducibility reasons. Investigation aspects are the possibility in identifying zones with uniform costs within the cost structure and the computation time.

III. GENERAL PROCESS OF THE FOR DETERMINATION BY THE PARTICLE SWARM OPTIMIZATION
In the following, the general process of the FOR determination (see Fig. 1) using a PSO-based aggregation method is introduced before considering costs for a flexibility provision by the FPUs in section IV.

A. DESCRIPTION OF FLEXIBILITY POTENTIALS AS PQ-POLYGON
The specific characteristics of a FPU are relevant for the FOR determination and must be communicated by the stakeholder to the corresponding system operator of a grid area (see I. in Fig. 1). The flexibility potentials of the FPU are nowadays limited by the regulatory framework of the technical guidelines. An example for the active and reactive power flexibility potential ( P min , P max , Q min , Q max ) of a throttled wind turbine based on the German technical guideline VDE-AR-4110 is shown in Fig. 4a) [44]. The flexibility potentials of the wind turbine are described as a polygon within the PQ-plane. In general the wind turbine manufacturers are guaranteeing a larger PQ-polygon than required [45], [46].
The flexibility potential of a FPU depends on the individual technology of the system element. The theoretical active and reactive power flexibility potential of a wind turbine with a DFIG is illustrated in Fig. 4b) [47], [48]. Compared to Fig. 4b) the technical guidelines (see Fig. 4a) lead to a loss of flexibility potentials for the system operator, which could be required and may be provided as ancillary service potential in specific system states. Differences between the technical requirements and the theoretical flexibility limits exist for each type of FPU. Fig. 5 shows six characteristic types for FPU PQ-polygons (cf. [2], [26], [30]) approximating their theoretical limits. These PQ-polygons were already used within case-studies in the context of FOR determination.
Even though the PQ-polygons in Fig. 5 represent the theoretical flexibility limits of the FPUs more than the requirements within the technical guidelines there are still deviations (cf. Fig. 4b) and type VI in Fig. 5). The first reason for this is the approximation of the unit circle for the apparent power through linearization. The second reason is the convexification of the PQ-polygons for a simple description of the flexibility potentials (see [2], [49]). In general, FPU PQ-polygons can be non-convex (see Fig. 4b) [49]. Therefore, methods for the FOR determination need to consider also non-convex PQ-polygons to maximize the utilization of the FPU flexibility potentials. To guarantee this for the methods developed within this article, the nonconvex PQ-polygon of Fig. 4b) is used for wind turbines with a DFIG in the case study in section V [47].
The general FPU PQ-polygons need to be adapted prior to the FOR determination according to the specific technology (e.g. P = 0 for compensation of Type I, minimum power supply at Type IV, on/off-states of a FPU) [2], [33], [47]. The flexibility potentials of the FPUs are influenced by the current operating point of the system element (cf. [35]). For example, a throttled operation is necessary for an increase of the active power supply of a wind turbine (see Fig. 4a). Further examples are the available primary power supply of a wind turbine or a PV-unit, the state of charge of a storage or incremental steps for the partial operation of a synchronous generator or load. For more realistic investigations, the wind forecast and the economic interests of the stakeholders can considered [35]. Nevertheless, a convex or non-convex FPU PQ-polygon results.
On load tap changing (OLTC) transformers represent another type of FPU, which can be used by the corresponding system operator for in-phase and/or quadrature voltage control [50]. OLTC transformers are not described by a PQ-polygon but instead by lower and upper boundaries and the incremental change of the voltage magnitude and angle per tap set.
In the context of coordinated DSO/DSO-interactions the flexibility potentials of a lower-level DSO system can be also described by a PQ-polygon representing a FOR (see I. in Fig. 1) [3], [16], [18], [51]. Thereby, the active and reactive power flexibilities of lower-level DSO system can be implemented analogously to PQ-polygons of FPUs by the higher-level DSO. Due to the non-linear system behavior of electric energy systems FOR PQ-polygons are in general nonconvex (see Fig. 3) [9], [23], [25], [41].

B. SAMPLING OF THE FOR EDGES BY THE PARTICLE SWARM OPTIMIZATION
The sampling of the FOR edges for a specific voltage V (see Fig. 1) is based on multiple solutions of the nonlinear OPF problem defined by the objective function in (1) due to variations of ϕ as well as the plus and minus signs [9], [23], [29]: By dividing the unit circle in ϕ steps l = 360 • / ϕ sampling processes results (see Fig. 6). Security constraints are represented by the minimum and maximum voltage limits, the maximum thermal current of the lines as well as the rated loading of the transformers.
The operational degrees of freedom are represented by the PQ-polygons of the FPUs. The population of the PSO is VOLUME 10, 2022 represented by a swarm consisting of n swarm particles [40], [52], [53]. The objective of the swarm is to find the minimum of the objective function through interactions of the swarm particles during their movement through the R 2k search space, where k is the number of buses. For each sampling process an individual PSO run is started. Thereby, each particle is described by a position and a velocity information. The position x of a particle represents the individual active and reactive power flexibility provision per bus by the FPUs: The velocity of a particle describes the change of the active and reactive power supplies of the FPUs for the next iteration step. At the beginning of the iterative solution process (t = 0) the swarm is initiated uniformly distributed within the limits of the FPUs PQ-polygons (see Fig. 6). Therefore, the FPU PQ-polygons are Delaunay triangulated by t triangles and described by the vertices vectors p P1 − p P3 and q P1 − q P3 : The cumulative sum a r,i of the triangle areas 1 to i related to the complete area of the PQ-polygon is given by (4).
Uniformly distributed, random numbers r = [r 1 , . . . , r j , . . . , r k ] in the interval of [0,1] are generated according to the swarm size k to identify a specific triangle for the initialization of the swarm particle positions. A specific triangle w is identified for each random number by: a r = [a r,0 = 0, a r,1 , . . . , a r,w−1 ≤ r i ≤ a r,w . . . , a r,t ] T (5) Two uniformly distributed random numbers r 1 and r 2 (r 1 ≤ r 2 ) within the interval of [0,1] (cf. [54]) are used to determine a random point within triangle w: During the convergence process of the PSO with a maximum number of t max iteration steps m = t max n power flow calculations based on the Newton Raphson algorithm considering the current swarm positions are performed to evaluate the objective function value. To avoid a swarm movement outside the FPU PQ-polygons a set-to-limit operator [55] based on the point-in-polygon test according to Jordan is used. Points outside the PQ-polygon are either set to a vertex or to the nearest point on the edges by orthogonal projection. Velocities, which lead to a movement outside the PQ-polygons are inverted for an improved convergence behavior. Particles that are not complying the technical constraints are punished and represent an invalid solution. For more details regarding the PSO algorithm (e.g. general algorithm, punishment function) used in this article see [23], [50]. The positions of the swarm particles and additional data from the power flow calculations (e.g. flexibility provision per FPU, bus voltages, active power losses) represent metadata.

C. AGGREGATION OF FPU FLEXIBILITIES BY THE MINKOWSKI SUM
In general, the performance and the quality of the results of optimization methods scale with the number of operational degrees of freedom and constraints [30], [40], [2], [56]. Furthermore, the consideration of multiple flexibilities contributing equally to the system (e.g. several FPUs at a single bus) can lead to a bad convergence behavior [49]. To avoid this, the complexity of the optimization problem is reduced by the aggregation (index: a) of the flexibility polygons F FPU,i of the n FPUs connected to a single bus (see Fig. 7a). A PQ-polygon F FPU i is represented by the coordinate vectors of the edges p FPU i and q FPU i : Prior to the aggregation, the flexibility polygons are moved to [0,0] (see F FPU 2 in Fig. 7a) according to the current operating points (index: op).
Therefore, the active and reactive power flexibilities of a FPU ( P min , P max , Q min , Q max ) are directly accessible:  For two convex flexibility polygons F FPU 1 and F FPU 2 the aggregated flexibility polygon F FPU,a is determined by the Minkowski sum (see Fig. 7b), which is the totality of all sums of the edges of the individual polygons [49]: Equation (9) is only applicable if one of the flexibility polygons is convex. If both flexibility polygons are nonconvex the area can be divided into convex sub-polygons (e.g. by Delaunay triangulation) [49], [57].

IV. MONETARIZATION OF THE FOR BASED ON A COST-OPTIMAL FLEXIBILITY DISAGGREGATION
The monetarization of the FOR is based on the principles of a local DSO and global TSO flexibility market [39], [58]. Within the local market the DSO has access to the flexibility potentials within the distribution system [37], [39], [58], [59].

A. ASSUMPTIONS FOR THE LOCAL DSO FLEXIBILITY MARKET
Today, the active power prize within the dayahead active power market is specified by a Market Clearing Price (MCP) settlement mechanism. The MCP corresponds to the maximum selling bid price. The provision of reactive power is mandatory according to the technical guidelines and is not remunerated [48]. Without an incentive regulation, stakeholders of generating FPUs are only interested in a maximum active power supply whereas stakeholders of loads are interested in an active and reactive power consumption according to demand. Stakeholder of storages are interested in low market prices for loading and high market prices for unloading.
In future, an availability and a utilization payment are necessary for the provision of ancillary services [48]. The availability payment (commodity prize in e/MW) arise from guaranteed accessible flexibility potentials. An example is to fund the throttled operation of wind turbines for positive secondary control power reserve. The utilization payment (service prize in e/MWh) results from the active and reactive power flexibility provision in a specific operating point. A local DSO flexibility market is characterized by individual, free bids of the stakeholders for commodity and separated service prices for the active and reactive power provision. Stakeholders of FPUs with short term preserved flexibility potentials can also attend at the flexibility market with the utilization payment mechanism. The availability payment mechanism is considered prior to the FOR determination to specify a pool of FPUs for guaranteed flexibility potentials. In general, aggregators are used to coordinate FPUs with small flexibility potentials but are neglected within this article.
At local DSO and global TSO flexibility markets the DSO defines its own flexibility demand and aggregates the remaining bids of the FPUs for the TSO [58]. Within this article, the DSO aggregates the flexibility potentials as monetarized FOR considering security constraints and costs for increased distribution grid losses [58]. The TSO has global access to the bids by the FPUs connected to the transmission systems as well as flexibility potentials within the FOR limits. For a certain flexibility demand the TSO specify the most economical adaptation of the IPF. The DSO is managing the TSO/DSO-coordination and has to guarantee this change during operational management [12], [58].
Considering commodity and utilization payment mechanisms two separated market processes are necessary. The commodity payment mechanism is less complex and can be adapted from the utilization payment mechanism, which is why it is not considered in the following.

B. SERVICE PRICES FOR THE ACTIVE AND REACTIVE POWER FLEXIBILITY PROVISION OF FPUs
A general concept for the determination of the reactive power flexibility service costs of a FPU is based on the expected payment function (EPF) in Fig. 8 [48], [60]. Three prize zones based on individual and free bids for linear cost factors (see I, II, III in Fig. 8) are used to consider the economic interests of a specific FPU stakeholder.
• Zone I: Fixed prize to meet the limits of the mandatory requirements from the technical guidelines.
• Zone II: Variable prize for the reactive power supply for example according to quadratic increasing converter losses (in Fig. 8: linear approximation) • Zone III: Variable prize for the reactive power supply according to a necessary active power reduction Different price zones within the PQ-polygon of FPU (e.g. zone I) are neglected for simplification reasons but an extension of the presented approach is possible. Because of the focus of the EPF to consider only reactive power flexibilities it cannot be applied without adaptations for the determination of a FPU service prize, considering active and reactive power flexibilities. Separated service prizes for the provision of active and reactive power flexibilities as combination of zone I and II are introduced.  Therefore, a linear monetarization of the active and reactive power flexibility provision P and Q for a duration d by the cost factors c s,P and c s,Q is defined: Depending on the specific FPU type (10) can be customized and different price zones of the FPU flexibility polygon may result. One example is the consideration of different costs for positive (+) and negative (-) active and reactive power changes: c s = P + ·c s,P,+ + P − ·c s,P,-+ Q + ·c s,Q,+ + Q − ·c s,Q,-d In the following, the service costs for the flexibility provision of the FPUs is based on (10) for an appropriate but simple description of the active and reactive power flexibility service costs.

C. SERVICE PRIZE DETERMINATION FOR EACH SWARM PARTICLE
The determination of the service prize c s,β of a single FPU connected to a bus β is already given by (10). In this case, the metadata of a swarm particle j for the flexibility provision per bus ( P β , Q β ) is used to determine the total service prize c s,tot,j for a specific vertical active and reactive power flow P vert and Q vert : In general, multiple FPUs with individual active and reactive power service costs and flexibility potentials are connected to a common bus β. Aggregated FPU PQ-polygons represent the summarized flexibility potentials (see section III C) to reduce the complexity of the optimization problem. The cost-structure of this aggregated FPU is unknown and (12) is not applicable. To enable the use of (12), the cost-structure of an aggregated FPU needs to be determined before the monetarization of the FOR. Within a premonetarization step the DSO disaggregates specific active and reactive power flexibility demands per bus ( P β , Q β ) cost-optimal to the individual FPUs to identify the coststructures of the aggregated FPUs. For the specification of P β and Q β a scatter within the limits of the aggregated FPU PQ-polygon is used as input for the MILP disaggregation described in the following.

D. DETERMINATION OF THE COST-STRUCTURE OF AGGREGATED FPUs BASED ON MILP DISAGGREGATION
The disaggregation of an active and reactive power flexibility supply to the single FPU within an aggregated FPU at bus β considering individual active and reactive power service costs represents a MILP in the general form of (13).
Each non-convex PQ-polygon is divided into a number of convex sub-polygons (e.g. Delaunay triangulation). Enhanced approaches can be used to reduce the number of convex sub-polygons and by this the complexity of the optimization problem. The variables P nc,µ , Q nc,µ describe a specific operating point within a FPU flexibility sub-polygon µ = 1, . . . , f . Analogously, the variables P c,η , Q c,η describe a specific operating point within the initial convex FPU flexibility polygon η = 1, . . . , g. The Boolean variable σ represents an either-or-condition to guarantee the utilization of a single sub-polygon per nonconvex FPU flexibility polygon. The variable λ is used to determine the absolute value of the service costs within the objective function. In (14) the inequality constraints of the MILP are presented: The matrices C nc , C c , B nc and the vector b c describe the limitation of the solution space by linear equations (m 1 · x 1 + m 2 · x 2 = b) based on the individual number of edges κ of 5408 VOLUME 10, 2022 each FPU flexibility polygon: The matrices I nc and I c are identity matrices of the size (f × f ) and (g × g), respectively. Line 3 to 4 in (14) belong to an if-then-condition to consider the σ corresponding constraints for non-convex FPU flexibility polygons which are divided into multiple convex sub-polygons. The matrix M max include the maximum active and reactive power values of the corresponding convex sub-polygon edges p µ and q µ : For the determination of the matrix M min the (17) is adapted for the minimum active and reactive power values of the corresponding convex sub-polygon edges p µ and q µ . Line 5 to 8 in (14) are used to determine the absolute value λ of the service costs (c s,P , c s,Q ) by using the cost matrices Co nc and Co c , e.g.: The equality constraints are given by: Lines 1 and 2 in (19) guarantee the compliance of the active and reactive power demand P β and Q β from the specific aggregated FPU flexibility polygon at bus β. I nc,µ and I c,η are identity matrices of size (2 × 2). The vector e assign the convex sub-polygons to the corresponding nonconvex, aggregated FPU flexibility polygon within an eitheror-condition. The minimum limits for the flexibility variables are given by: x min = min( p 1 ), min( q 1 ) · · · min( p f ), min( q f ) × min( p 1 ), min( q 1 ) · · · min( p g ), min( q g ) × 0 · · · 0 0 · · · 0 0 · · · 0) T (20) The maximum limits of the flexibility variables are given in (21).
× c s,P,1 max( p 1 ), c s,Q,1 max( q 1 ) · · · × c s,P,g max( p g ), c s,Q,g max( q g ) T (21) The total service prize c s,tot,j for each particle j can be determined according to (12) based on the presented MILP disaggregation method. An example for the determination of the cost-structure of an aggregated FPU within the premonetarization step is presented within the case study in section V.

E. DETERMINATION OF THE FOR COST STRUCTURE
An overlay of the monetarized t max particle swarms lead to overlapping points within the cost structure of the FOR (see Fig. 3a) and an identification of cost zones is not possible. Overlapping points represent swarm particles with similar values for P vert and Q vert but different flexibility provisions per bus x = [ P 1 , . . . , P k , Q 1 , . . . , Q k ]. For the TSO only the particle with the minimum summarized service prize of the FPUs is relevant for a cost-optimal flexibility provision. The costs for the active power losses P loss,j of particle j and a duration d are determined by (22) and added to c s,tot,j : c loss,j = P loss,j − P loss,0 c loss d ∀ P loss,j − P loss,0 > 0 (22) Within (22) only increased grid losses compared to the initial losses P loss,0 are monetarized for the TSO. The idea behind this is that a reduction of grid losses is in the economic interest of the DSO.
For the determination of the FOR cost-structure, the swarm particles within the whole iteration process of the PSO are assigned to a scatter within the limits of the FOR PQ-polygon. The scatter is described by the active and reactive power vectors p sc and q sc (see Fig. 9). For the determination of the cost structure c sc the Euclidean distance between the active and reactive power interconnection power flow of a swarm particle (P vert,j , Q vert,j ) and the scatter points (P sc,y , Q sc,y ) is identified (see Fig. 9a): For the scatter point y with the minimum distance to the position of a swarm particle j the cost value c sc,y is updated by the particle costs c s,tot,j , if: c sc,y c sc,y > c s,tot,j = c s,tot,j Based on (24) the particles with the minimum costs for the corresponding scatter points are identified (see Fig. 9b). The monetarization of the FOR is completed by the determination of the cost structure.

V. CASE STUDY
An adaptation of the Cigré medium voltage system is used for the case study [42]. The original grid was extended by aggregated low voltage grids and a variety of different FPUs. Each type of FPU is working at the same operating point and has the same flexibility potentials adapted from [2]. In scenario 1, the FPU PQ-polygons correspond to Fig. 10. In scenario 2 the installed power of the FPUs is halved. The high voltage bus is used as slack with a specified voltage of V slack = 110 kV e j0 • . Information on the lines and the transformers are given in Table 1 and 2. The voltage change per tap-set of the OLTC HV/MV-transformer at the higher-voltage bus is 0.25% related to the corresponding rated voltage. The maximum and minimum tap-set is limited to ±10 steps. Service costs for the OLTC transformer are neglected within this article. In general, the service prize of a OLTC transformer can be estimated by the equivalent lifetime loss per tap-set [17]. The full data set, including the FPU flexibility polygons, is available at [43].    The number of samples for the FOR determination is 45 with a corresponding sampling angle of ϕ = 8 • in (1). The swarm size n and the maximum iteration step t max are set to 200, which leads to 1.8 million power flow calculations. For a detailed parameterization of the PSO see [23]. The sampling process is fully parallelized on 45 workers with usual computation power in MathWorks MATLAB.
The prize c loss in (22) is set to 50 e/MWh according to typical German market prices. The service cost factor (see Table 3 and 4) for the active power flexibility provision The objective of the case study is the application and evaluation of the process for the determination of the FOR cost structure, presented in section IV. The investigations are divided into two steps. First, the cost structure of the aggregated FPU PQ-polygon is premonetarized by the MILP presented in section IV D. Exemplarily, the five FPUs at bus 2 and the FPU PQ-polygons of scenario 1 are used. Second, in both investigation scenarios, the monetarized FOR at the vertical system interface is determined based on the methods described in sections IV C and IV E.

A. PREMONETARIZATION OF AGGREGATED FPU PQ-POLYGON
The composition of the PQ-polygon of the aggregated FPU at bus 2 is presented in Fig. 11a). Starting from the load (see I in Fig. 10) the flexibility potentials of the other FPUs are added by the Minkowski sum. The PQ-polygon of the aggregated FPU is non-convex due to the DFIG flexibility potentials. The MILP disaggregation (see section III D) is applied for a 100 × 100 scatter in the limits of the aggregated FPU PQ-polygon for the determination of the cost structure (see Fig. 11b). Further points are specified by the edges of each subpolygon combination during the Minkowski process in Fig. 11a).
The linear service cost factors c s,P and c s,Q of the FPUs at bus 2 leads to a two-dimensional piecewise linear function for the cost structure of the aggregated FPU. The information of the cost structure can be used for a cost-optimal distribution of an active and reactive power flexibility demand at bus 2 on the individual FPUs. The computation time is 0.9 s.
For the premonetarization of the three FPUs at bus 5 the computation time was 0.5 s and for the busses with two FPUs about 0.2 s. The computation time depends on the complexity of the MILP and the number of MILP runs, which depends on the resolution of the meshgrid. With a higher number of FPUs within the aggregated FPU the computation time for the premonetarization will increase. To avoid this, the FPUs can be sorted and aggregated regarding equal economic interests of the stakeholders prior to the premonetarization. The computation will be also increased for solving a mixed-integer quadratic programming problem in the case of quadratic cost functions for the FPUs active and reactive power service prizes. These aspects become important at more realistic grid scenarios with a variety of different FPUs at a single bus.

B. DETERMINATION OF THE FOR COST STRUCTURE IN SCENARIO 1
Within the heat map of Fig. 12a) the density of the swarm particles during the PSO-based sampling process and the resulting FOR are presented. The edges of the FOR are sampled more detailed then the center of the FOR, which results from the convergence behavior of the PSO (cf. Fig. 6). The FOR is limited by the technical constraints of the grid (see Fig. 12b) and not by the flexibility limits of the FPUs.
The results of the FOR including the service-costs of the FPUs and the costs for increased grid losses within the distribution grid are shown in Fig. 13a) and b) based on a 100 × 100 scatter. Within the FOR the cost zones are blurred, which results from smooth cost gradients. Exceptions are the areas near the maximum active power supply and reactive power consumption of the lower-level system. Considering the density of the swarm in Fig. 12a), these areas are sampled by fewer swarm particles, which result in less metadata.
From this it can be derived that for areas with a low swarm particle density in consequence of a challenging system state (here: high voltages, high transformer loading) the cost structure becomes less uniform. Another reason is that the flexibility potentials of the FPUs are larger than the maximum power VOLUME 10, 2022 transfer capability of the HV/MV-transformer. Thereby, more constellations of the FPUs to guarantee a specific IPF are possible. The monetarization of the FOR is only based on metadata because the reduction of the service costs is not an objective within the sampling process of the PSO. This results in regions with varying service costs (see Fig. 13a) in the usually uniform cost-structure. The non-uniform transition of the cost structure in the middle of the FOR results from the lower density of the swarm particles in this area.  Within the case-study only higher active power losses compared to the initial system state are monetarized in (22). Thereby, the dark blue area in the bottom, middle and top of the FOR in Fig. 13b) result. The grid losses are increased especially in case of high active power transfers at the vertical system interconnection in case of high load or supply. The transition of the service costs within the complete coststructure of the FOR (see Fig. 14) is similar to the results in Fig. 13a). The impact of the active power losses is only significant in the area with a high active power transfer to/from the lower-level system. In Fig. 15a) and b) the cost structure of the FOR is presented in more detail. The point with the lowest service prizes is represented by the initial system state without any flexibility provision by the FPUs. Especially for particles located in the middle of the FOR the flexibility provision of the FPUs is not representing the minimum service costs. In general, the higher-level system operator can assume the results for the cost structure as estimation for the real costs to be excepted for a specific vertical active and reactive power flow. The computation time for the FOR sampling is 56 s and for the determination of the cost-structure 6 s. In scenario 2 the density of the swarm particles within the FOR is higher than in scenario 1 (see Fig. 16a). The top and bottom edges of the FOR are limited, analogously to scenario 1, by the lower and upper voltage limits (see Fig. 16b). The left and the right side of the FOR are limited by the flexibility limits of the FPU PQ-polygons.  Within the cost-structure of the FOR in Fig. 17 clear cost zones can be identified. This also applies for regions with a low density of the swarm particles like the maximum active power consumption (see Fig. 16a). The variations of the cost gradients are more significant compared to scenario 1. This can be identified by the width of the individual color ranges.
The computation time for the FOR sampling in scenario 2 is 58 s, which can be explained by numerical performance variations. The computation time for the determination of the cost-structure is 8 s. The increased computation time results from the higher number of swarm particles complying the technical constraints within the FOR compared to scenario 1. Within scenario 1 only 1.2 million swarm particles are included in the FOR compared to 1.6 million in scenario 2.

VI. DISCUSSION AND CONCLUSION
The presented approach for the monetarization of the FOR leads to comprehensible results within the case-study of the Cigré medium voltage system and the additional computation time for processing the metadata is small. In both investigation scenarios, cost-structures with mostly clear cost zones are obtained. In general, the cost structure is blurred in regions with a low sampling density. Significantly more definable cost zones can be achieved by combining multiple costs to a common zone (e.g. 40 e−50 e, cf. [14]). The computation time for the premonetarization of the aggregated FPUs and the processing of the metadata will increase for larger systems with a variety of buses and FPUs. Investigations on larger systems are necessary to evaluate the performance of the algorithms and the quality of the results in more detail. Nevertheless, the presented approach is suitable for an application within the dayahead and intraday operational management.
The cost-structures are representing an estimation of the service costs that can be assumed by the higher-level system operator for a specific adaptation of the IPF and the corresponding flexibility provision of the FPUs. The costoptimal distribution of a flexibility demand (e.g. ancillary service provision) from the higher-level system operator to the FPUs within the lower-level system is guaranteed within the distribution step of the hierarchical grid control strategy (see III. in Fig. 1). This can result in lower service costs for the higher-level system operator. Nevertheless, the plausibility of the cost-structure can be investigated in more detail by comparing the costs for specific, representative operating points with the results of a cost-optimal provision of the corresponding IPFs. Based on the results an improved sampling of the FOR regarding monetarization aspects can be developed to reduce the gap.

VII. SUMMARY AND FUTURE RESEARCH ASPECTS
This article continues research regarding the aggregation of ancillary service flexibility potentials within lower system levels at the vertical system interconnections in the context of hierarchical multi-(voltage-)level grid control strategies and especially TSO/DSO as well as DSO/DSO-cooperation. The main contribution is the monetarization of the active and reactive interconnection power flows (IPF) within the feasible operation region (FOR) by metadata from the particle swarm optimization (PSO). Advantages of the presented approach are the simple integration within the previous aggregation process without an extended sampling process. Furthermore, the possibility to consider a variety of different flexibility providing units (FPU) with individual economic interests for commodity and service prizes and even non-convex active and reactive power flexibility polygons is advantageous. Therefore, the aggregated flexibility potentials of multiple FPUs connected to a common bus are premonetarized by a cost-optimal flexibility disaggregation to the individual FPUs within a mixed-integer linear programming problem (MILP) prior to the FOR determination. Thereby, the cost-structure of the aggregated FPU can be determined in any level of detail. The metadata generated by the PSO during the sampling of the FOR and the resulting cost-structure of the FPUs are used to specify the service costs for a specific IPF and to monetarize the FOR. By this, an enhanced approach for a techno-economic TSO/DSO-coordination based on a local DSO and a global TSO flexibility market results compared to current literature.
The promising results of the case study are the basis for further extensions of the presented process and the application of other aggregation methods like random search and other metaheuristics. The potential extensions can be divided into three categories.
The first category is coping with the implementation of further functionalities to the developed approach and investigations regarding new research aspects: • Specification of more realistic investigation scenarios (e.g. the shape of the flexibility polygons is influenced by partial load operation and on/off-states of the FPUs besides the regulating of the FPUs) • Consideration of technology specific FPU cost factors • Development of an enhanced local DSO and global TSO flexibility market mechanism considering the FOR • Extension of the aggregation method to multiple vertical system interconnections • Consideration of (n-1) scenarios The second category adapts FOR related research aspects from literature to implement them within the process presented in section III and IV. Examples are the consideration of uncertainties due to forecast deviations, time constants for the FPU flexibility provision or the use of voltage controlling transformers as additional flexibility potential. VOLUME 10, 2022 The third category deals with the further development and implementation of the monetarized FOR within the cascading process of a hierarchical multi-level grid control strategies. Therefore, the monetarized FOR at the DSO/DSO system interface has to be integrated within the operational management of the higher-level DSO. Beside a specification of an individual flexibility demand by the higher-level DSO the next step is the aggregation of the flexibilities for the TSO at the TSO/DSO system interface. Therefore, the methods and processes presented in this article serve as a promising basis.