Flexibility Services Management Under Uncertainties for Power Distribution Systems: Stochastic Scheduling and Predictive Real-Time Dispatch

Over the last years, the role of the distribution system operator (DSO) has largely expanded. This is necessitated by the increased penetration of intermittent energy resources at the distribution level, as well as the new, more complex interactions with the transmission system operator (TSO). As such, to properly manage its system and to have an effective joint cooperation with the TSO, the DSO is required to procure and carefully manage flexibility services from distributed energy resources (DER). This paper introduces a thorough framework on optimal operational planning (day-ahead scheduling) and operational management (real-time dispatch) of active distribution systems under uncertainties, to avoid line congestions and voltage limit violations, and efficiently balance the distribution system. A two-stage stochastic programming model based on weighted scenarios is proposed to optimize the multi-period optimal power flow day-ahead scheduling, i.e., scheduled power flows at the TSO-DSO interface and reserved DER flexibility services. Subsequently, the operational management, realized with a predictive real-time dispatch model based on a constantly updated rolling horizon, aims to efficiently activate the available flexibility services to minimize deviations from the committed schedule. Different sources of flexibility are considered, with their respective response times also taken into account at real-time dispatch. The proposed framework is applied on two distribution systems and investigates the DSO’s level of risk exposure while minimizing its total cost (reservation and activation expenses). The results indicate that a less conservative approach at planning stage, despite the potential risk exposure, can lead to significant reduction in total expenses.


I. INTRODUCTION
The ever-increasing integration of renewable DER due to the deregulation of the energy market creates more challenges for the DSO, as high levels of uncertainty are introduced into the electrical power system, especially at the distribution level. Apart from the traditional problem of voltage control in passive distribution systems, the issue of congestion management has also attracted significant attention. In adhering with the future smart grid vision [1], active distribution systems are expected to host very high shares of renewable technologies, which have been noted to cause substantial thermal issues if left uncontrolled. It comes to no surprise that recent research works have proposed approaches to combat congestions [2]- [6]. On the other hand, DER are also excellent sources of flexibility for managing the electrical systems. The DSO can utilize DER flexibility to optimize its system's voltage profile, manage local congestions [7], and generally optimize the operation of smart distribution networks through ANM schemes [8].
However, the uncoordinated activation of DER flexibility by the DSO is likely to impact the system's balance and imply extra costs to the TSO, as explained in [7]. A prerequisite for the beneficial use of DER flexibility is the joint cooperation with the TSO, a topic that, recently, is widely discussed by many studies [9]- [11], with several coordination schemes having been proposed [12]. For the purposes of this paper, the TSO-DSO coordination is realized under the scheme of shared balancing responsibility model, according to the conceptual framework of [12]. Under this coordination scheme, the DSO organizes a local flexibility services market and is responsible for congestion management and balancing the distribution grid, while staying committed to a predefined schedule in the interconnection with the TSO.
This paper aims to propose a comprehensive framework of flexibility services management for the DSO. Within this framework, the DSO firstly carries out the DAS of the necessary DER flexibility for congestion management and balancing the local system, while optimizes an import/export hourly schedule at the T-D to timely notify it to the TSO. Subsequently, the distribution system is operated by proper RTD of the reserved DER flexibility. The overall objective is to minimize the deviations from the committed import/export schedule at the T-D with minimum flexibility services costs (reserve and activation cost), while meeting all technical limits.
In the past, the DAS was modelled as a multi-period OPF problem, assuming negligible forecast error, and thus ignoring the inherent uncertainties of RES and loads. Articles [13]- [15] are representative of such a DAS approach presenting deterministic models that anticipate the periods when network technical constraints are violated and propose volt/var control and/or up/down active power regulation for congestion management during these periods. In [13], the DAS is modelled as an active-reactive OPF that considers the four-quadrant operation capability of the BSS to maximize wind energy yield. A dynamic OPF formulation with energy storage and flexible demand in an ANM framework is proposed in [14]. An integer genetic algorithm is employed to optimally solve the DAS problem and minimize operational costs in [15].
In recent years, research has focused on new models that integrate in DAS the uncertainties of RES and load. The authors of [16] propose a deterministic, dynamic OPF model that considers the operational uncertainties of wind generation and demand during the planning stage. This model schedules the import/export at the T-D and commits reserves for the worst-case uncertainty scenario. Such an approach leads to more expensive reserves than an approach that would consider more scenarios of uncertainty along with their probability of occurrence. In [17], a two-stage stochastic multi-period OPF is used to schedule the generation of the dispatchable units given the available demand response flexibility, while considering the wind power uncertainty under multiple weighted scenarios. The authors of [18] introduce a two-stage robust OPF model for DSOs willing to utilize flexibility services from local DER to solve congestion/voltage problems considering the worst-case scenario of uncertainty. Even though the approach of [18] concludes in more reliable solutions, it also leads to more expensive reserves and with total operational costs similar to those of a deterministic scheduling with OPF.
The DAS is by itself insufficient to properly manage the distribution system, and thus, RTD is additionally employed to achieve improved operation. Several dispatch techniques have been proposed, from single-period OPF models to more advanced techniques of predictive control based on multiperiod OPF models with a rolling horizon. RTD with predictive control is based on the principles of MPC, a model that predicts the future evolution of a process to optimize a set of control variables over a constantly updated rolling horizon.
Several studies have employed the MPC approach to manage the operation of power systems [19]- [23]. In [19], a MPC algorithm is proposed to solve the economic dispatch problem in the presence of RES in electrical systems, while the authors of [20] propose a centralized MPC strategy to maximize the wind energy yield in an isolated electrical system by compensating for wind and load variability. Studies [19]- [20] both highlight the potential of MPC in electrical systems, but overlook important technical aspects, such as voltage and power flow limits, without having formulated the power flow equations. A combined day-ahead and intra-day hierarchical optimization is introduced in the works of [21]- [22] to maximize RES exploitation, while minimizing variations between DAS and RTD. A MPC-based dispatch approach addresses the uncertainties of electric vehicle charging patterns in [23].
However, none of the proposed models has considered the response time of DER offering flexibility services, which is proposed as a new flexibility index by [24]. The response time is also presented by [25] as a physical requirement of flexibility that should be taken into account by the operators. Both [24] and [25], highlight the significance of integrating the response time in the modelling of flexibility services. To the authors' knowledge, this is the first paper to take into account the response time of flexibility services in the MPC algorithm to solve the RTD problem, with respect to DSO utilizing flexibility services from DER.
The contribution of this paper is threefold: a) It proposes a novel framework for managing flexibility services that optimizes both the DAS of the reserved flexibility as well as their activation during RTD procedure. The proposed framework additionally investigates the DSO's VOLUME 8, 2020 level of risk exposure that minimizes the total expenditures of flexibility services (DAS and RTD costs).
b) It proposes a novel stochastic DAS approach based on a two-stage stochastic programming formulation. The proposed S-DAS schedules the power flows at the T-D and minimizes the reserved flexibility services by minimizing the cost of their expected activation.
c) A predictive control strategy is employed for the RTD procedure. The proposed P-RTD not only considers the modelling of the ramping limits and the dynamic operation of DER, but also proposes a model for the efficient handling of the different response times of system's flexibility sources.
The remainder of this paper is structured as follows: Section II presents the proposed framework for flexibility services management. Section III formulates the proposed S-DAS. Section IV describes the proposed P-RTD optimization procedure. Section V presents the examined case studies and analyses the obtained results. Section VI concludes the paper and proposes future research work.

II. PROPOSED FRAMEWORK FOR FLEXIBILITY SERVICES MANAGEMENT
This paper proposes a framework for optimal operational planning (DAS stage) and optimal operational management (RTD stage) of active distribution systems by means of congestion management and system balancing, while maximizing the active power injection from RES. With this approach, the DSO exploits the DER flexibility services in the most effective way in order to stay committed to the hourly import/export schedule at the T-D with minimum deviations. The expected value of these deviations is minimized at the DAS stage by anticipating the needed flexibility services, whereas the actual value of these deviations is minimized during RTD by activating the needed portions of flexibility services from those reserved during DAS.

A. MAIN ASSUMPTIONS
For the sake of clarity, the main assumptions adopted in this paper are summarized as follows: a) The FSP included in this study are: i) DERA (aggregators of flexible loads and DER operating at LV), ii) BSS (operating at MV), and iii) RES (operating at MV). Their models are described in Section III.
b) The proposed framework is applied under normal operating conditions, where the scenarios considered are solely related to the intermittent nature of RES and load. Extreme conditions, such as WG tripping, which require contingency reserves, are not within the scope of this paper. c) Wind power generation depends, among other factors, on various parameters of weather conditions. In this paper, for simplicity, the authors have considered that wind power generation depends only on wind speed and thus, is calculated based on the well-known wind turbine S-curve P w (v) [13]. Additionally, without loss of generality, wind speed profile is considered the same across all buses of the distribution system.
d) It is supposed that an online local flexibility market platform has been developed. This platform is intended to facilitate the trading between DSO and the FSP.
e) The FSP are able to submit their bids in the local flexibility market platform until 10:00 am. After the local market is cleared, the DSO timely notifies the scheduled power flows at the T-D to the TSO, not later than 12.00 at noon. The TSO, in turn, clears the ancillary services day-ahead market under its responsibility, having the predefined schedules at the T-D.

B. FLEXIBILITY SERVICE PROVIDERS
A flexibility service, as defined in [26], is a power regulation (upward or downward) sustained at a given time-period for a predefined duration from a specific location within the network. The flexibility service, as a product, is technically characterized by four indices: a) power capability, b) energy capacity, c) ramping limit and d) response time [24].
In what follows, flexibility services are split in two main categories according to the direction of active power regulation. Thus, the term FSD refers to demand reduction (or generation increase) and vice versa for the term FSU.

1) DER AGGREGATORS
A DERA is an intermediary entity that aggregates the demand of end-customers' flexible loads and the generation of residential DERs (e.g., roof-top photovoltaics, other small-scale RES) to provide both FSD and FSU to the DSO upon request. The DERA provides its services at the MV level of the distribution system. The DSO, in turn, offers proper remuneration, according to the DERA's bid (energy-price pair), depending on the reserved and the activated volumes of FSD and/or FSU. Therefore, the DERA is remunerated at a ''commitment price'' (in e/MW-h) for keeping the needed FSD and FSU reserves, and at an ''activation price'' (in e/MWh) for the finally delivered FSU and FSD.
FSD offered by DERA can be achieved by reducing customers' aggregated consumption, or by increasing the power injected by local DER (e.g., discharging of the scheduled stored energy for backup). FSU, in accordance, is derived by increasing customers' aggregated consumption (e.g., through price signals or demand response programs) or by reducing residential RES aggregated production (e.g., generation curtailment). The way in which the DERA achieve FSU and FSD is not the focus of this paper.

2) BATTERY STORAGE SYSTEMS
BSS can inherently provide bi-directional flexibility services, both FSD and FSU, when discharging and charging, respectively. It is assumed that BSS are privately-owned, and an agreement is in place stating that the DSO may reserve partial of full capacity of the BSS at the day-ahead stage. The BSS are then dispatched by the DSO to ''absorb'' any power mismatches from the day-ahead committed schedule. The portion of BSS's capacity that is scheduled to be utilized is determined during DAS. In case the capacity of the BSS is partially reserved, it is assumed that its owner uses the remaining capacity in a way that does not affect the corresponding bus's net load profile. Here, again a twofold remuneration scheme is in place for the services provided by BSS owners. The BSS owner is remunerated by the DSO at a ''commitment price'' (in e/MW-day) for its capacity availability and at an ''activation price'' (in e/MWh) when charging/discharging active power is activated.

3) RENEWABLE ENERGY SOURCES
Privately-owned large-scale RES connected at MV of the distribution system are subject to curtailment providing a service of FSU. RES curtailment is the DSO's last resort, as it is usually activated in case of forecasted congestion or when energy mismatches between DAS and RTD are more expensive than RES tariffs. The conditions under which a DSO can force curtailment are determined in a power purchase agreement contract between the DSO and the RES owner. Curtailment policies, as well as the compensation paid to the RES owners differ per country [27]. In this paper, the cost of RES curtailment is set equal to a fixed price (in e/MWh), supposing a feed-in-tariff contract.
Large-scale RES should also be able to provide reactive power support according to grid code. The reactive power requirements, as specified in different grid codes, are summarized in [28]. For the purposes of this study, the reactive power capability ranges ±33% of the active power injected by RES, i.e., power factor ranges from 0.95 leading to 0.95 lagging. This is a producer's obligation and thus the DSO uses reactive power compensation free of remuneration.

1) DAY-AHEAD SCHEDULING (1 ST STAGE)
Formulated as a multi-period OPF, in the general approach, DAS comprises the decisions to be made one day in advance given the day-ahead forecasting. In this paper, a two-stage stochastic programming model for the multi-period OPF, hereafter called S-DAS, is proposed to optimize the first stage variables, while ensuring feasible control actions throughout the second stage. To simulate the uncertain realization of the second stage, multiple scenarios are considered, based on the day-ahead forecast, weighted according to their probability of occurrence.
The first stage variables, pertaining to here-and-now decisions that ensure feasible solution across all scenarios considered, are the following: 1) scheduled import at T-D, i.e., the HV/MV SS, 2) reserved FSD and FSU by DERA, 3) reserved BSS power capability. The vector of here-and-now control variables in time-interval t is defined as follows: The RTD stage is the validation process of the DAS. In this stage, the DSO transitions to the RTD, given the reserves and the schedule of DAS. In general, RTD is an iterative procedure that determines the control actions to be implemented at each time-interval of the operation based on the decisions made in the DAS stage and given the updated state of the system. In this paper, a P-RTD is proposed, formulated as multi-period OPF to take advantage of a look-ahead rolling VOLUME 8, 2020 horizon, based on the principles of MPC. The modelling of response time of the available flexibility services is also introduced. The P-RTD utilizes the constantly refreshed shortterm forecasting to determine the real-time control actions, anticipating optimal future control actions of the rolling horizon.
The P-RTD, thoroughly analysed in Section IV, optimizes the vectors of the second stage variables u P−RTD t+h|t and x P−RTD t+h|t at the h−th time-interval of the t−th rolling horizon, given the system's current state, i.e., for h = 0 of the previous rolling horizon, u P−RTD t|t−1 and x P−RTD t|t−1 .

III. S-DAS FORMULATION
In this section, firstly the uncertainties of RES and load are modelled by creating different scenarios. Subsequently, the models of the FSP, the interaction with the TSO and the network constraints are described. The objective function and the overall S-DAS formulation are also given.

A. DAY-AHEAD FORECASTING UNCERTAINTIES MODELLING
This subsection explains the process of scenario creation for handling the uncertainties of the day-ahead point forecasts of load and RES generation. These scenarios represent scenarios of the distribution system's net load, since they are derived from combinations of different uncertain states of load and RES generation. Their values are used as input parameters for the S-DAS.

1) PROBABILISTIC FORECASTS WITH PREDICTION INTERVALS
In decision-making, probabilistic forecasting provides a better view of the potential future outcomes than deterministic forecasting. The term deterministic forecasting refers to methods for producing point forecasts, which represent the conditional expectations of the future outcomes without considering future errors, whereas the term probabilistic forecasting is used when prediction error and its probability density function are considered. A prediction interval is a type of probabilistic forecast directly linked to past errors due to point forecasts, which do not consider the error distributions [29]. According to [30], a prediction interval with a coverage rate α% is defined as a range of potential values for the future observations y t+h , computed based on past observations, such that the future observations y t+h will fall in this interval with probability α, where α ∈ [0, 1]. For example, for a prediction interval with a coverage rate of 90%, there is 90% probability for the future observations to belong in this prediction interval [30]. Prediction intervals are usually based on the root mean square error of the deterministic forecasting model, which provides an estimate of the standard deviation of the forecast error [31]. In this paper, the upper and lower bounds of the prediction interval are constructed by (1) assuming that the forecast error of the h − th next time-interval is normally distributed with zero mean and standard deviation σ h , i.e., ε h ∼ N 0, σ 2 h [32].
whereŷ t+h denotes the point forecast issued at time-interval t for the h−th next time-interval and z is the critical value of the standard normal distribution that gives the coverage rate of the prediction interval. For example, z = 3 gives a prediction interval with a coverage rate of 99.7% [31]. Do note that the standard deviation of the forecast error depends on the lead time of the prediction.
In what follows, the probabilistic forecasts of the future outcomes are realized based on a prediction interval with a certain coverage rate. This prediction interval, in turn, is used to create multiple representative states of the future outcomes, each one with its probability of occurrence. These states are created by dividing the selected part of the normal distribution of the forecast error (i.e., the area within the prediction interval selected) into distinct intervals, according to the methodology followed in [33] and [34].
Wind speed and load forecast errors are both assumed to follow the normal distribution, a common approach employed by many researchers to model such forecast errors [34]- [36]. Modifications on the probability distribution function of the forecast errors could be made accordingly to better approximate the reality; however, this is out of the scope of this paper. Nonetheless, the forecast error modelling (for both wind speed and load) that is detailed hereafter is generic; it is applicable for any distribution function desired, should one wish to obtain even more representative results.

2) WIND SPEED FORECASTING UNCERTAINTIES
The uncertainties of wind speed deterministic forecasting are strongly related to the forecast error of the numerical weather prediction tools. This error is assumed to be normally distributed with zero mean and standard deviation σ w [35].
To construct multiple states of wind speed, firstly, the normal distribution of the wind speed forecast error is divided into distinct intervals. An example of nine intervals is presented in Table 1. The corresponding value of wind speed forecast error for each state (ε w v ) is determined taking into account the worst-case per state in accordance with Table 1.
Following this process, a set of states for the wind speed forecast error is created, each state with its value and its probability of occurrence: {ε w v , π v }.

3) LOAD FORECASTING UNCERTAINTIES
The DSO should also predict the day-ahead demand at each of the system's buses. Assuming that the load forecast error follows the normal distribution, the same approach with wind speed probabilistic forecasting is applied. Likewise, the example of Table 1 is applied for the states of load forecast error. The probability of each state is calculated with (3a), while the sum of all probabilities of the selected states is equal to the coverage rate a L of the considered prediction interval, as shown by (3b). A set of states for the load forecast error is created, each state with its value and its probability of occurrence:

4) SCENARIOS CREATION
For the creation of scenarios, each state of wind speed forecast error is combined with each state of load forecast error, i.e., {ε w s , ε d s } = {ε w v , ε d l } s , using a scenario tree model as in [37]. Assuming that the individual states of wind speed and load forecast error are independent, the probability of occurrence of scenario s is calculated according to (4a). The load forecasting curve under scenario s is calculated by (4b) based on the day-ahead point forecastP frc d,t . The wind power forecasting curve under scenario s is given by (4c), as a function of the wind speed day-ahead point forecastv frc t by using the S-curve transformation of the wind turbine of WG w, i.e., P w (v). Equation (4d) gives the sum of all probabilities of the created scenarios. The value of β represents the probability the upcoming observations of load and wind speed to be covered by the created scenarios. Hence, β is the coverage rate of the selected scenarios.
The participation of DERA in the day-ahead flexibility market is described by (5a)-(5e). Equations (5a) and (5b) denote the maximum energy of FSU and FSD the DSO can reserve over a day. The power capability of FSU and FSD are formulated by (5c) and (5d), accordingly. Equation (5e) gives the cost of reserves that should be remunerated.
The operation of BSS under all scenarios is modelled by (8a)-(8f). The BSS energy balance at each time-interval should satisfy the dynamic mathematical formulas (8a)-(8b). The upper limit of power that can be absorbed by or injected to the system is described by (8c)-(8d). Equation (8e) models the upper and lower technical limits of the BSS. Equation (8f) gives the activation cost of BSS charging/discharging.
The penalties for upward and/or downward deviations are set by the TSO and are related to the real-time system price and the cost of extra energy portions the TSO must dispatch to balance the system under its responsibility.

D. AC POWER FLOW MODEL
The AC power flow equations, given by (11a)-(11c), are applied ∀i, t, s and ∀j ∈ Ni : The node power balance is expressed by (12a)-(12b): Equations (13a)-(13c) describe the system's technical constraints. Equation (13a) sets the bus voltage limits. Equations (13b) and (13c) set the capacity limits of network's power flows and HV/MV SS, respectively. The objective function of S-DAS, given by (14a), consists of two terms, each one related with the stage the decisions refer to. The first term is (14b) and the second term is the expected value of (14c), both over the 24-hour scheduling The general mathematical formulation of the S-DAS is given in (15) as follows: subject to The proposed S-DAS problem is a NLP optimization problem and is solved using the state-of-the-art CONOPT4 commercial solver of GAMS [38]. In S-DAS, the duration of time-intervals is considered equal to one hour (d t = 1), i.e., 24 time-intervals.
The S-DAS problem is inherently a large-scale optimization problem. Its computational complexity primarily depends on: a) the number of network buses, b) the multiple periods of the day-ahead planning and c) the number of the scenarios that model the uncertainties of wind speed and load. For instance, to optimize the day-ahead planning (24 hourly periods) for a 33-bus distribution system considering 9 states of wind speed and 9 states of load (81 scenarios of net load), the time elapsed is approximately 34 minutes. All tests in this paper were performed on a PC with an Intel Core i7-8700 CPU at 3.20 GHz and 8 GB of RAM. When the number of buses increases, the tractability degrades in a nonlinear fashion.
The computational performance could be improved by employing convex relaxations of the AC power flow equations, such as the semi-definite programming relaxation [39], the second order cone relaxation [40]- [41], and the quadratic convex relaxation [42].
In this paper, the authors' main interest is to propose a novel and effective S-DAS model to solve the problem and validate the extracted results, rather than improve the execution performance. Given the novelty of the proposed approach, it is paramount to ensure its reliability, as is done in this paper, before developing further add-ons.

IV. P-RTD FORMULATION
In this section, firstly a model to handle the response time of the FSP is introduced. Subsequently, the algorithm of the proposed P-RTD procedure is outlined based on the principles of MPC. Moreover, the overall P-RTD formulation is given.

A. HANDLING RESPONSE TIME
The time delay between receiving the activation signal sent by the DSO and achieving the requested output is defined as the response time of the flexibility service. The response time is a critical parameter during RTD as it plays significant role in the optimal flexibility sources management and the efficient activation of the less expensive flexibility services on time. Response time varies from a few seconds to few minutes and depends on the technologies used to implement the flexibility services, as analyzed in [24] and [25]. The handling of the response time during P-RTD can be depicted by Fig. 2. Whenever the response time is equal to the timeslot of the dispatch procedure, the vector of the current control actions u t n |t n cannot be immediately dispatched and is replaced by the vector of the control actions u t n+1 |t n−1 of the second timeslot of the previous rolling horizon. In case the response time is equal to twice the timeslot of the dispatch procedure, the vector of the current control actions u t n |t n is replaced by the control actions u t n+2 |t n−1 of the third timeslot of the previous rolling horizon, and so on. For technologies with response time quite faster than the timeslot of the dispatch procedure, it is assumed that the vector of the current control actions u t n |t n is applied at the beginning of the current timeslot.
BSS are characterized by a fast response time, up to some seconds, and can quickly compensate mismatches attributed to short-term forecasting errors. RES curtailment can be activated with response time up to some seconds, but the response time may be set to longer durations (e.g., some minutes) depending on system operator rules to avoid unwilling ramp-up production, e.g., due to wind gusts. Here, RES maximum permitted production is announced 15 minutes in advance. As far as the response time of DERA is concerned, the activation signal is usually sent from 5 to 30 minutes in advance depending on the energy technologies used and their contractual agreement with the DSO. In this paper, a response time of 15 minutes is assumed, which is a prerequisite for many DSOs when procuring flexibility services [43].

B. PROPOSED P-RTD BASED ON MPC PRINCIPLES
In general, MPC allows for the optimization of current timeslots, while always taking future time-slots into account, by constantly re-predicting future events based on the actual flow of real events. This is achieved by continuously optimizing the system's operation under a finite horizon but implementing the control actions only on the current operating condition, thus performing a brand-new optimization for every single operating condition, under an ever-changing rolling time horizon.
In the proposed P-RTD, the DSO solves a look-ahead multi-period OPF problem for each time-interval. The optimization depends on the current value of the system's dynamic state variables, i.e., BSS state of charge at the end of the previous time-interval and the DERA's set-point of the previous time-interval due to ramping limits. The procedure of the proposed P-RTD is implemented as follows: a) at time-interval t the multi-period OPF problem is solved over a predefined rolling horizon with H time-intervals. d) the optimization process is repeated from step a) for the next time-interval t = t+1 having new short-term forecasting values for the next rolling horizon.
Based on the above procedure, the control actions of BSS charging/discharging and DERA's FSU and FSD are constantly updated to minimize the objective function of (18a) over the finite rolling horizon by considering: a) the most recent forecasted data, b) the control actions of the flexibility sources with greater response time, and c) the current state of the system.

C. P-RTD OVERALL FORMULATION
Equations (6), (8)- (13) are used in this stage by excluding the subscript s (where it exists) and by replacing the subscript t with t + h|t that indicates the variable of the h − th timeinterval of the t − th rolling horizon. Equation (9a) is reformulated with (16a)-(16b), where P Inj w,t+1|t−1 is the set-point (maximum permitted injection) of RES w for h = 1 that was optimized in the previous rolling horizon (t − 1).
In P-RTD, (14c) is modified to formulate the objective function of the t − th rolling horizon, as follows: The general mathematical formulation of the multi-period OPF problem of the t − th rolling horizon is given by (18). The iterative procedure of P-RTD has been coded in MATLAB, while the optimization model (formulated as a NLP problem) is solved using the state-of-the-art CONOPT4 commercial solver of GAMS [38]. In P-RTD, the duration of time-intervals is considered less than one hour with d t ∈ (0, 1], e.g., if the duration is 15 min, d t = 0.25 and the number of rolling horizons over a day is 96.
In contrast with the computational complexity of S-DAS, the P-RTD of each time-interval is executed in less than 10 seconds. This is mainly due to the absence of the scenarios of uncertainty in P-RTD.

V. CASE STUDIES & NUMERICAL RESULTS
This section is dedicated to investigating the DSO's level of risk exposure by applying the proposed flexibility services management framework. At the very beginning of this section, the methodology that evaluates the risk exposure is introduced. Subsequently, the proposed flexibility services management framework is implemented in two case studies. The first is based on a 4-bus distribution system and is presented as a proof-of-concept validation. The second is based on a 33-bus distribution system and is presented for further investigation.

A. RISK EXPOSURE EVALUATION METHODOLOGY 1) CASES OF RISK EXPOSURE
Four different cases (presented in Table 2) are considered in order to investigate the DSO's level of risk exposure. Each case of risk exposure is linked to a certain prediction interval of load (α L ) and a certain prediction interval of wind speed (α V ), the coverage rates of which are considered equal, as presented in the second column of Table 2. The third column of Table 2 gives the prediction interval bounds, for both load and wind speed, following the notation of (1) and assuming that both errors are normally distributed with zero mean. The prediction intervals of load and wind speed are selected based on the states of Table 1, which are shown in the fourth column of Table 2. Pertaining to the last column, the rate of risk exposure is 1 − β, with β ∈ [0, 1], where β is the coverage rate of the selected scenarios as defined in (4d) of Section III.A.4. The selection of the cases of Table 2 is based on the modelling of Section III.A.
Case A is the most conservative among the selected cases with a rate of risk exposure (1−β) of only 0.6%, which means virtually zero risk exposure for the DSO. To create the scenarios of Case A, all the states of load forecast error and wind speed forecast error that are shown in Table 1 of Section III.A are included, the combination of which concludes in 81 scenarios, each with its own probability of occurrence.
On the other hand, Case B, Case C and Case D demonstrate the DSO's level of risk exposure given that the higher the rate of risk exposure (1 − β) is, the less the reserves are, and thus the riskier the RTD is. Therefore, Case D is the least conservative among the selected cases, because it has a rate of risk exposure (1 − β) of 53.3%. Transitioning from Case A to Case B, the states that lie on the tails of the forecast error's normal distribution, i.e., those with smaller probability (1 st and 9 th states) are excluded. Combining the 7 remaining states of load forecast error with those 7 remaining states of wind speed forecast error, the 49 scenarios of Case B are created. Their probabilities are normalized to neglect the probabilities of the rejected scenarios in S-DAS, following the theory of the truncated normal distribution applied in [44]. Following the same procedure and excluding more states, Case C and Case D comprise 25 and 9 scenarios, respectively.

2) RISK EXPOSURE METRIC
The metric additional risk exposure, denoted by ARE X →Y , is introduced to measure the risk exposure, when transitioning from Case X to Case Y, given that Case X is more conservative (has lower rate of risk exposure) than Case Y. As shown in (19), ARE is the ratio of the increase in RTD cost over the reduction of DAS cost, when transitioning from Case X to Case Y. In fact, ARE is a ratio that measures the increase in RTD cost as a result of reducing the DAS cost by an additional monetary unit (e.g., e).
where RTD X and DAS X denote the RTD cost and the DAS cost of Case X, respectively. The values of ARE can be interpreted as follows. If ARE X →Y has a value less than unit, the transition from Case X to Case Y is beneficial for the DSO, since the DAS cost reduction is higher than the RTD cost increase. Inversely, if the value of ARE X →Y is greater than unit, it is no worth transitioning from Case X to Case Y. Figs. 4 and 5 depict the probabilistic forecasts of network's total load and wind generation, respectively, with their corresponding prediction intervals (shadowed areas) in accordance with the cases of Table 2. The point forecasts and the actual measurements are also shown to illustrate the actual deviations from the day-ahead point forecasts.
With regard to the remuneration schemes described in Section II.B, Table 3 shows the prices concerning FSP   participation in this test system. The penalties for deviating from the import schedule at the T-D are given in Fig. 6 (the penalties of upward deviation are considered equal to those of downward deviation).
For the day-ahead forecasting, the standard deviation for the load forecast error and that for the wind speed forecast error are considered equal to 10% and 15%, respectively, of the point forecasts [16]. It is usual the wind speed forecast error to be greater than that of the load forecast [45].
As far as the short-term forecasting is concerned, the predictor of the P-RTD procedure is assumed to have an error with standard deviation that varies over time, from 3% (1 h ahead) to 15% (24 h ahead) for the wind speed, and from 1% (1 h ahead) to 10% (24 h ahead) for the load forecast (expressed as percentage of the point forecasts).

2) S-DAS RESULTS
The results of S-DAS, presented in Table 4, show that the more conservative the DSO is, the higher the reserve costs become. The total DAS cost of Case A (e545.9) with rate of risk exposure of only 0.6% is more than twice as expensive compared to the total DAS cost of Case D (e245.0) with rate of risk exposure of 53.3%. Most of the cost that is avoided by adopting riskier approaches is attributed to the significant reduction in the level of DERA reserves, while the reservation costs of the remaining flexibility options, i.e., the BSS, remain relatively at the same level, as Table 4 shows. This makes DERA the additional option in guaranteeing system balance, as well as the source for most of the DSO's dayahead expenses. Fig. 7 shows how the total energy of FSD is allocated across the day, making clear that when considering Case C and Case D for S-DAS, the reserved FSD energy is much lower than that of considering Case A and Case B. The narrower the coverage rates of the prediction intervals of load and wind speed are, the less the expected deviations from the point forecasts are supposed to be, and thus lower FSD reserves are needed.

3) P-RTD RESULTS
Using the reserves scheduled with the S-DAS, the DSO applies the proposed P-RTD to stay committed to the hourly import schedule at the T-D, while exploiting the available flexibility services with the least possible expenses. The P-RTD is executed with time-resolution of 15 minutes and is tested for different rolling horizons to highlight the importance of look-ahead optimization in RTD.    Fig. 8, the first set of bars, labeled with ''no P-RTD'', refers to an ordinary RTD strategy that does not consider look-ahead horizon. The set of bars ''no P-RTD'' is used as benchmark for comparison. The rest sets of bars refer to the length of the selected rolling horizon of the P-RTD, e.g., the label ''4h P-RTD'' denotes P-RTD of a 4 h rolling horizon. Fig. 8 shows that in all cases of risk exposure, the proposed P-RTD has better performance (lower RTD cost) than ''no P-RTD''. The longer the rolling horizon of the P-RTD is, the lower the RTD cost is.
Depending on the length of the rolling horizon, the P-RTD process differs in accuracy (the extension of the rolling horizon increases the forecast errors) and thus different costs are incurred. For instance, it is observed from Fig. 8 that the RTD cost of Case A is reduced from e649.5 (''no P-RTD'') to e555.6 (''1h P-RTD''), i.e., a reduction of 14.5% from ''no P-RTD''; when extending the rolling horizon by 1 hour the RTD cost of Case A is reduced from e555.6 (''1h P-RTD'') to e527.1 (''2h P-RTD''), i.e., an additional reduction of 5.1%; and when extending the rolling horizon by 2 hours (''4h P-RTD'') the additional reduction from ''2h P-RTD'' to ''4h P-RTD'' is about 1%. While the extension of the lookahead horizon generally reduces the costs incurred, the additional reduction quickly diminishes. Extending the rolling horizon's length makes sense until a certain point. After a certain threshold is reached, it appears that the marginal gains are minimal. In fact, if the P-RTD approach is also coupled with a very broad look-ahead horizon, it may lead to additional costs, due to the fact that the more extended is the horizon the higher is the variance of the forecast error. Thus, apart to how much value the DSO gains from considering prediction intervals with higher coverage rate for S-DAS, there is also a case to be made on how much value the DSO gains from extending the rolling horizon of P-RTD.
There is merit to take a closer look at Fig. 8 and analyze the reduction of RTD cost for different cases of risk exposure, while keeping the length of the rolling horizon to a fixed value. For example, focusing on the set of bars ''4h P-RTD'', the RTD cost of Case A is reduced from e649.5 (''no P-RTD'') to e520.1, i.e., a reduction of 19.9% from that of ''no P-RTD''; the RTD cost of Case B (equal to that of Case A) is also reduced by 19.9%, while the RTD cost of Case C and Case D are reduced by 17.5% and 10.9%, respectively, from those of ''no P-RTD''. This shows that the performance of P-RTD, not only depends on the rolling horizon selected, but also depends on the flexibility services reserved, highlighting the dynamic nature of the P-RTD.  Figs. 9, 10 and 11 refer to P-RTD of 4 h rolling horizon (''4h P-RTD'') and present results of RTD when considering Case C for the S-DAS. The results of Figs. 9, 10 and 11 are compared with the results obtained if a dispatch strategy without predictive control (''no P-RTD'') is implemented. Fig. 9 illustrates the import at the T-D and compares the curves of ''no P-RTD'' and ''4h P-RTD'' with that scheduled during DAS. This comparison determines the deviations of actual active power flows from the committed ones at the T-D. Fig. 10 depicts the FSD activation that the DERA provides, while Fig. 11 presents the operation of the BSS by depicting its state of charge. By simultaneously analyzing Figs. 9, 10 and 11, one can notice that P-RTD ends up in lower deviations during the time window 16:00-21:00 (between 64-th and 84-th time-interval), when T-D deviation penalties (Fig. 6) are much higher, than during early morning hours.

4) DAS & RTD OVERALL DISCUSSION
The results extracted for the illustrative example (Sections V.B.2 and V.B.3) signify the vital need for careful S-DAS. In the pursuit of reducing the reserve costs, the adopted risk may be very high, and may lead to hefty intraday balancing costs that are not offset by the savings of the DAS stage, regardless of how sophisticated the utilized real-time prediction tool is. Fig. 12 presents the total cost breakdown in RTD cost and DAS cost for the four cases of risk exposure defined in Table 2 of Section V.A. In Fig. 12, the RTD cost refers to P-RTD simulation with a 4 h rolling horizon. Table 5 gives the values of ARE metric (defined in Section V.A.2). As an example, the ARE for the transition from Case B to Case C (denoted as B→C in Table 5) is computed as follows: The ARE A→B = 0.00 (Table 5) shows that the transition from Case A to Case B is beneficial for the DSO, because no extra RTD cost is incurred as a result of the additional gains achieved through DAS cost reduction by e40.8. Increasing once more the rate of risk exposure (1 − β) and transitioning from Case B to Case C, ARE B→C = 0.30. This indicates that the B→C transition reduces the DAS cost more than three times compared to the RTD cost increase, as is also shown in Table 5. Consequently, the transition B→C is also beneficial for the DSO.
The ARE C→D = 2.03 (Table 5) shows that the transition from Case C to Case D is not beneficial for the DSO. Indeed, although the DAS cost is very low in Case D (Fig. 12), the DSO has severely jeopardized system balance, and as a result, the RTD cost of Case D is much higher than that of Case C (Fig. 12). The risk undertaken in Case D (much higher RTD cost) is not properly compensated by the DAS cost reduction. Keeping less reserves in Case D causes much   higher T-D deviations (from the schedule), which are accompanied by stiff penalties (those of Fig. 6), as well as more expensive flexibility activation.

C. CASE STUDY OF 33-BUS DISTRIBUTION SYSTEM 1) CHARACTERISTICS OF THE DISTRIBUTION SYSTEM
The 33-bus distribution system of [46], modified according to Fig. 13, is also studied. The basic network data are directly sourced from [46]. This network is radial, with nominal voltage of 12.66 kV, and is supplied by a 10 MVA capacity HV/MV substation. The maximum line capacity is set at 8.7 MVA and the peak load is considered equal to 4.37 MVA. The network contains controllable loads by DERA, and also has multiple installed DER: WG, PV stations, and grid-scale BSS. The WG total installed capacity is 3.2 MW and the PV station's installed capacity is 0.6 MWp. Two grid-scale BSS of 1.3 MW total power capability and 2.5 MWh capacity (90% round-trip efficiency) are also part of the system. Two DERA are competing in this case study, by regulating the net load of different buses: DERA-1 controls    accordingly. The prices remunerated to the FSP that participate in this case study are given in Table 6.

2) RESULTS
Fig. 14 shows the breakdown of total cost for the four cases of risk exposure presented in Table 2 of Section V.A, for the 33-bus distribution system. Fig. 14 shows that Case C results in an overall cost (DAS and RTD cost) of e971.72, which is the lowest value among the selected cases. Even though the rate of risk exposure is higher when deciding to solve S-DAS with Case C, the DSO reduces its overall cost by 7.3% in comparison with that of Case A (e1,048.62), which is the most conservative case. Once again, there is a threshold of risk exposure, over which the DSO's total expenses begin to increase. Indeed, as is shown in Table 7, the values of ARE metric are less than unit for the transitions A→ B and B→C, while ARE is greater than unit for the transition C→D.
Interesting findings can be drawn when multiple DERA are participating in the local flexibility market. Fig. 15 shows the reserved FSU and FSD versus those activated by the two DERA participating in the 33-bus distribution system, i.e., DERA-1 and DERA-2. The results of Fig. 15 refer to the reserves of FSU and FSD considering Case C for the S-DAS. The activation of FSU and FSD considers a 4 h rolling horizon for P-RTD. As Fig. 15 shows, the DSO activates more the FSD of DERA-1 (4.13 MWh) than that of DERA-2 (0.45 MWh). However, with a closer inspection of Fig. 15, the degree the DSO exploits the reserved FSD of DERA-2 (63% of the reserved FSD, which is 0.72 MWh) is higher than that of DERA-1 (48% of the reserved FSD). This is reasonable, because the activation price the DERA-2 has bidden is cheaper than that of DERA-1. Of course, DERA-1 FSD total energy is higher, as the energy bid is also much higher than that of DERA-2. However, when it comes to FSU, even though DERA-2 activation price is still lower, the FSU reserves are zero owing to the higher commitment price (e25/MW-h) compared to that of DERA-1 (e20/MW-h), and thus no FSU could be activated by DERA-2.

D. OVERALL DISCUSSION
There is a clear trade-off between the reserved flexibility services and the real-time expenses. In other words, sacrificing a part of the reserves of Case A, the DSO spends less in total (DAS plus RTD expenses). What is more important is to investigate the maximum level of risk exposure, i.e., the portion of the reserves that could be sacrificed leading to lower total expenses. This varies among different test systems and depends on the cost of reserves, the cost of activation and the penalties of deviating from the scheduled power flows at the T-D. Further investigation is needed in this direction, in order to be able to elaborate more the results of different test systems for different range of net load scenarios to draw more specific conclusions for the flexibility services management under uncertainties.

VI. CONCLUSION AND FUTURE WORK
This paper introduces a flexibility services management framework consisting of two optimization stages: 1) a stochastic scheduling approach to efficiently handle both dayahead uncertainties and intra-day fluctuations, and 2) a realtime dispatch based on predictive control used as validation of the DAS. The work addresses contemporary issues for both optimization stages that are frequently disregarded. The general insufficiency of the deterministic day-ahead scheduling is demonstrated due to the always present forecasting errors. Even so, the inclusion of stochastic programming to model the planning stage can prove inadequate, leading to over-or underestimations of the planning decision due to the high variability of renewable technologies. In addition, the dire need for a dynamic and highly adaptable approach for real-time dispatch applications is demonstrated by employing the MPC principles. The constantly updated short-term predictions on the output of wind farms are such that require constant monitoring with a continuous stream of updated setpoints timely sent to the various rapidly responding technologies (e.g., BSS), which are a powerful tool in the DSO's arsenal.
With respect to the recorded numerical results, it has been shown that when performing stochastic programming, there is a limit to the level of risk exposure, i.e., what range of scenarios should be considered, when compared to the marginal gains, especially if they are not high enough to offset the extra RTD cost. While additional reserved flexibility is certainly not undesirable, approaches that are overly conservative are generally avoided in the industry sector, as they usually lead to more expensive solutions. The limits of risk exposure with respect to marginal gains are illustrated by the various values of ARE metric that is introduced to measure the DSO's level of risk exposure. It is also shown that when performing optimization for real-time dispatch based on a look-ahead technique, there is a saturation to the effectiveness of extending the time horizon.
The handling the response time of the FSP is also a critical aspect, which is illustrated in this paper. It is the first time that a model is formulated to effectively activate the less expensive flexibility services on time, while considering a look-ahead technique. The quantification of the benefits from handling the response time is also valuable, but not within the scope of this paper; however it is strongly recommended as a future work.
Moreover, an issue not within the scope of the paper but with significant research interest is the tractability of the proposed approach. While proven effective, its reliable application in much larger distribution systems is challenging. As such, the proper approximations must be devised to expand the applicability of the proposed approach for larger problem sizes. Furthermore, an issue that has been ignored by past research works is the case of the TSO making emergency requests for flexibility from the DSO. The approach that the DSO would have to adopt is to include these new unplanned (time-wise and quantity-wise) parameters in its schedule, and it is an issue that has not been properly addressed in the scientific literature. It is in that light that the authors plan to significantly extend their work by not only incorporating even more modern aspects of active distribution networks, but also doing so in a scalable fashion.