Multimarket Trading Strategy of a Hydropower Producer Considering Active-Time Duration: A Distributional Regression Approach

This article presents a new approach for finding the optimal multimarket trading strategy of cascaded hydropower plants (HPPs) in the sequential electricity markets. These markets are day-ahead energy market, the market for frequency containment reserve in normal mode (FCR-N), and manual frequency restoration reserve markets for both energy production and capacity reserve. The active-time duration (ATD) of an mFRR energy offer is an important required parameter and it is uncertain at the time of day-ahead offer-function submission. Hence, we suggest a distributional regression approach for ATD modeling in an optimal multimarket setup. Also, a modified machine learning approach is proposed to generate price scenarios for the mFRR energy market taking into account uncertain ATD parameters. To illustrate our proposed approach, various numerical experiments are performed. Our numerical results show how proper modeling of ATD parameters can lead to a more realistic multimarket offer-function for cascaded HPPs. Furthermore, the results show how the inclusion of FCR-N and mFRR capacity markets change the optimal day-ahead offer-function.

Probability of occurrence of each price scenario in second stage. π ω Probability of occurrence of each price scenario in third stage. k m Scaling coefficient for value of stored water calculation.

I. INTRODUCTION
A. Motivation T HE optimal schedules in the forward electricity markets such as day-ahead and intraday markets do not often ensure the strict operational security of the power system in real time. The transmission system operators (TSOs) keep the power network within its secure limits to ensure the security of supply at all times. The green transition toward a climate-neutral, secure, and integrated energy system, as the main vision of the most European TSOs [1], requires the new electricity-market designs for managing arising transmission-network related challenges [2], [3]. Flexibility service providers play an important role for a successful green transition and the hydropower plants (HPPs) have considerable benefit in providing flexibility services (such as ramp-rate and frequency-regulation services). However, the HPPs' technical constraints and their cascaded topology introduce some challenges in the planning and operation stages (e.g., information privacy [4] and operational constraints [5]). One of these challenges is the value of stored water which makes the HPPs along a river dependent on time and watercourse. The stored water in each reservoir has an opportunity cost (water value), making its reserve capacity procurement dependent on the available water in the downstream reservoirs.
In the European electricity market, there are several newly announced markets to address the challenges confronted by the TSOs to balance the system [3]. The frequency during normal operation is kept within ± 0.1 Hz using primary control. Below 49.9 Hz, the disturbance reserve is activated. The procurement markets consist of the frequency containment reserve for normal condition (FCR-N) and disturbance condition in European terminology. Like tertiary market in U.S., there is manual frequency restoration reserve (mFRR) energy market to retrieve the frequency to the nominal value after stabilizing the frequency due to the sudden change in the grid. Also, there will be a new capacity market for mFRR aimed to ensure TSOs that they have enough capacity available to be activated in the mFRR energy market during the operation day [6], [7]. In some countries, such as Germany and Austria, the mFRR capacity and energy markets are procured jointly once per day through a pay-as-bid remuneration pricing rule [8] and it is shown in [9] that splitting these two markets and procuring the mFRR energy market based on the marginal pricing reduces the balancing costs even in the presence of strategic players. Recently, the EU guideline on electricity balancing announced a common framework to initiate splitting the procurement of mFRR capacity and energy reserves [10]. In such multimarket setup, the optimal coordinated offer-function of HPPs is a complex task that requires a proper decision-making framework.
B. Literature Review 1) Optimal Offer-Function in the Multimarket Setup: Optimal offer-function calculation in the multimarket setting has been addressed in the previous literature with different assumptions and perspectives [11], [12], [13], [14], [15], [16], [17]. In [11], a detailed coordinated offering strategy in the day-ahead energy and balancing markets has been presented in which the electricity price and dispatched volumes are uncertain at the time offer-function submission. A qualitative discussion on offerfunction calculation in the hydrodominated multimarket setup has been performed in [12] using different mathematical optimization approaches. In [14], the extra benefit from coordinated offering in the day-ahead and balancing energy markets over sequential offering in hydrodominated power system has been estimated between 0.8% and 2.6%. However, in all the previous literature, it is assumed that the whole accepted balancing energy offer will be activated, which is not realistic as the TSOs will activate those offers for a short period of time which in turn influences the expected revenue of units.
Power producers benefit from trading in the reserve markets from the following two aspects: 1) they will receive compensations for being reserved regardless of whether the reserved capacity is activated, and 2) they can earn money proportional to the increase or decrease of their production when the reserves are activated. This is a profitable situation for HPPs as one of the flexibility sources from the TSO point of view. This participation can release some pressure from HPP operators who are struggling to recover their capital costs as well as their original profitability in the forward and spot markets partly due to the eroding peak/off-peak spreads [18]. Also, it is shown in [19] that allocation of 10% of transmission line capacity for exchange of reserves between bidding zones in the Nordic power market would result in a benefit of 290 k€ per day. However, naive consideration of the mFRR energy offer ATD leads to the unrealistic results and provides the possibility of energy arbitrage between day-ahead energy market and mFRR energy market (e.g., [20]), which is illegal from regulatory frameworks [21].
2) Statistical Methods and Scenario Generation Techniques: The electricity price information is often publicly available, and there are different approaches to generate price scenarios for the associated multistage stochastic optimization [22] models. The accuracy of the uncertainty representation depends on the number of scenarios and available information about the probability density function (PDF) representing the uncertain parameter. However, the active-time duration (ATD) of mFRR energy offers is uncertain at time of the offer-function calculation. By ATD, we mean the length of time in which TSO activates the procured mFRR energy Offers. This is an important factor that changes the optimal offer-function and expected revenue from the operators' point of view because their remuneration in the balancing energy markets is dependent on the ATD. Several factors, such as the net load imbalance, intermittent nature of the renewable energy generation, and forecast errors influence the ATD.
Generalized additive model for location, scale and shape (GAMLSS) is a well-known distributional regression approach, which has been recently developed and used in different applications [23], [24]. In GAMLSS, one of a variety of density functions can be selected in which the mean, variance and higher statistical moments can be specified dynamically as conditional upon exogenous driving variables. The authors in [25] and [26] have employed GAMLSS to model dynamically varying distribution of prices. Thus, we have used GAMLSS for modeling the uncertain ATD of mFRR balancing energy offers. Using the GAMLSS model, we show that the assumption of normal distribution for the ATD as the random parameter as proposed in the previous works (e.g., [27]) is not valid.
Furthermore, there is a need to generate electricity energy price scenarios for calculating the optimal multimarket offerfunctions. These scenarios need to be derived from a suitable price-forecasting model. Existing electricity forecasting methods can be divided into the following three categories, namely, physical methods [28], statistical methods (such as seasonal autoregressive integrated moving average (SARIMA) [29]), and machine learning methods [30], [31], [32]. Recently, with the development of deep learning, new electricity price prediction models have been continuously proposed [22], [30], [33]. A long short-term memory (LSTM) network is a special kind of deep learning models that is capable of learning long-term dependencies. LSTM deep neural networks have a good performance in handling nonlinear and complex dependencies in data. This method is suitable for cases in which the sequence of data is important. That is why its applicability for the electricity price forecasting in which different types of seasonality exist are recommended (see the literature review in [34]). The day-ahead electricity market price forecast based on RNN and LSTM has been studied in [35]. Also, in [30], the authors discuss the impacts of Nordic market coupling on the electricity price prediction through employing different feature selection techniques on the LSTM networks.

C. Technical Contributions
Our article contributes to the relevant literature as follows. First: On the modeling aspects, FCR-N and mFRR capacity markets have been considered as two revenue alternatives for HPP operators. The mFRR capacity market constraints limit the offer volume in the mFRR energy market as the operators should be ready to be dispatched with full volume submitted earlier in the capacity market. However, the FCR-N market does not bind the real-time energy market as it is only a capacity market (its energy remuneration is negligible). We have addressed the interaction of these two capacity markets in our proposed formulation.
Second: we have proposed a three-stage stochastic optimization model that finds the optimal day-ahead energy market offer-function of the HPPs along a river participating in four different market platforms, including the proposed mFRR capacity markets in the first contribution.
Third: Despite previous literature which have considered full ATD of mFRR energy offers or like [27] which employed the uniform and normal distribution to model the reserve deployment and wind power generation, we propose a distributional regression approach using the GAMLSS model with dynamic moments to find the best fit for the uncertain ATD. We then investigate the impact of our GAMLSS-based ATD model on accuracy of the suggested stochastic program.
Fourth: A forecasting method based on the LSTM deep learning architecture and GAMLSS model has been developed to generate scenarios for the third stage of our proposed stochastic program. This method utilizes the existing data for balancing energy markets from the second stage scenarios to update the third stage electricity price scenario forecasts. Those second stage scenarios have been extracted from historical data. In fact, we have updated the network state of the trained LSTM model with available data instead of predicted ones to achieve more accurate scenarios.
In a nutshell, in the problem formulation of an HPP operator in the hydrodominated electricity market participating in the new balancing markets, It is important to observe the interaction of these new market setups and possible offering strategies of the operator facing them (first contribution). To explain this interaction, we have proposed a new three-stage optimization problem that requires information about different scenarios at each stage (second contribution). Thus, we have proposed a modified machine learning-based scenario-generation technique to achieve better forecasts for the scenarios (fourth contribution). Meanwhile, we have noticed that the ATD of balancing energy offers in the real time, activated by the TSO, is an important factor in the decision-making framework that affects the expected revenue and requires appropriate modeling. Thus, we have proposed a GAMLSS-based approach to find the best distributional model for the ATDs (third contribution).
The rest of this article is organized as follows. Section II elaborates on the proposed multistage offer-function optimization problem and its constraints. Section III describes the GAMLSS model to find the best PDF fit for the uncertain ATD and proposes the LSTM-based price forecasting technique. Section IV provides the numerical results based on a Nordic hydropower case study. Finally, Section V concludes this article.

II. PROPOSED STOCHASTIC OPTIMIZATION FOR MULTIMARKET OFFER-FUNCTION CALCULATION
We assume a portfolio of reservoirs and their connected generation units owned by a company, which their operator offers into the day-ahead energy, FCR-N, mFRR capacity, and energy markets. The optimal trading strategy of this operator can be obtained through a three-stage stochastic program set out in (1).
In this formulation, we consider offering to the day-ahead energy market in which an HPP operator submits an aggregated offer curve for every hour of the following operation day. Here, we assume that offer prices are fixed as parameters while the volumes dispatched based on the revealed prices are decision variables. This assumption makes the optimization problem linear and solvable by off-the-self solvers like CPLEX and Gurobi. For a given hour t and seg = 1, . . ., N seg , the supply curve is defined by the prices ρ seg,t in which ρ seg−1,t ≤ ρ seg,t and ρ 1,t = 0, ρ N seg+1 ,t = +∞. The offer volumes, x spot seg,t ≥ 0, are assumed to be the accumulated volumes at a specific offer price. It should be noted that submitting offer-functions to FCR-N, mFRR capacity, and energy markets are out of the scope of this work and we consider that all of the offers in these markets will be accepted by TSO. In the following, formulation, we have t ∈ T , s ∈ S, and ω ∈ Ω subject to : The objective function (1a) maximizes the expected revenue from selling the electricity to different markets considering the value of stored water. First, two terms are the revenues from participation in the upward and downward mFRR capacity market. Third term is revenue out of FCR-N market. The fourth term is the revenue out of day-ahead energy market. Fifth and sixth terms are the revenue from mFRR energy market in the upward and downward trends and the last term is the value of stored water in the downstream reservoirs. Due to the large number in the value of stored water calculation, we have deducted a portion of the initial reservoir level as a constant value to make the revenue streams from different markets more comparable. Offer-function constraints and its nondecreasing features are stated in (1b)-(1c) and (1d)-(1e). The TSO requirements for the minimum volume of mFRR capacity are imposed in (1f)-(1g). Constraints (1h) and (1i) show the minimum and maximum capacity available for the FCR-N market considering the volume allocated to the mFRR capacity requirements and dispatched day-ahead energy market. Balancing energy market constraints, which limit the offer volume to the mFRR energy market with respect to the procured FCR-N and dispatched day-ahead energy market volumes are expressed in (1j)-(1k). The hydrological constraint on the reservoir level is presented in (1n) and energy balance is set in (1p). It is important to note that the impacts of ATD on the revenue stream and energy balance are considered through T up ω,s,t and T down ω,s,t . The piece-wise linear mathematical representation of hydropower output with respect to the discharged level is demonstrated in (1p). Constraints (1q)-(1s) enforce the operational limits on the variables based on the available data for each power plant. It should be noted that no decision is made between the time which day-ahead energy prices and FCR-N prices are realized. Hence, there is no need to create additional branching in the scenario tree between these events. We instead model pairs of day-ahead energy and FCR-N prices in the same scenario set s ∈ S. In this formulation, we have three stages elaborated as follows.
1) First stage: At this stage, the decisions for the offer volume in the day-ahead energy and mFRR capacity markets are decided. Furthermore, the HPP operator needs to allocate some part of its available capacity to the FCR-N market given the expected price of this market. It should be noted that there is no available information about the price of the mFRR capacity market as it is supposed to be launched by Q4, 2022 [6]. So, we have considered the mFRR capacity prices as the weighted average of mFRR energy market (λ dir t = k dir s∈S ω∈Ω π ω π s λ dir ω,s,t ), which dir represents the direction of balancing market which can be up or down.
2) Second stage: At this stage, the price information for day-ahead energy and FCR-N markets are revealed. Based on the offers submitted in the previous stage and the new information, the dispatched volume for each unit is determined. Furthermore, HPP operator needs to decide on the mFRR energy market offers before gate closure time [3]. 3) Third stage: At this stage, the price information for mFRR energy market is revealed and the HPP operator needs to decide how the operational variables (discharge, spillage, reservoir level, and electricity generation) should be determined. The graphical illustration of the abovementioned stages is depicted in Fig. 1.

A. LSTM-Based Scenario Generation Technique
In this section, we want to generate electricity price scenarios for the second and third stage of our formulation in (1). For the second stage, we have extracted one year of historical data and then reduced the number of scenarios using the fast-forward algorithm in GAMS/SCENRED [36]. For the third stage scenario generation, we have used a novel recurrent neural networks (RNNs) method. RNNs have been known as highly effective systems for a variety of sequence learning problems. An RNN can be considered as a generalization of feed-forward neural networks to sequential data, e.g., text data or time-series data [37]. To alleviate such training-related complications, Hochreiter and Schmidhuber [38] introduced the LSTM unit or block. The principal idea behind an LSTM block is a modifiable memory cell that can store information over long periods. The state of this cell can be altered by various gating units, allowing the network to learn what information is relevant and "unlearn" information that has become obsolete with respect to the purpose at hand. Naturally, the features described above facilitate the implementation of LSTM-based architectures in a wide array of problems involving sequential data.
In this article, we modify the scenarios generation techniques through LSTM architecture by updating the network state with available values instead of using the previously forecasted ones. This can happen if we have access to the actual values of times steps between predictions. In our case, after selecting the day-ahead energy market scenarios, their corresponding mFRR energy market prices are available (solid black lines in third stage of Fig. 1). But we need to generate required number of scenarios for third-stage (N 3 rd req ) given the available data for one scenario. The flowchart of this approach is depicted in Fig. 2. In this approach, X train is the one year of historical data for the mFRR energy market in each direction (upward and downward) [39] andX n,train is the normalized one-hour lagged version of the training data. X test is the test data set of mFRR energy balancing market, which is available based on the selected day-ahead energy market in the second stage, which is used to improve the forecasting. In the normalize section, we standardize the training data to have zero mean and unit variance. The training stage of the LSTM network is based on a sequence-to-sequence regression LSTM network, in which the responses are the training sequences that their values shifted by one time step (one-hour lag-generator in Fig. 2). That is, the LSTM network learns to predict the value of the next time step at each time step of the input sequence. The numerical analysis on this method is explained in Sections IV-B and IV-E. Table I shows the settings of the LSTM training network.

B. GAMLSS-Based PDF Estimation
In GAMLSS, the exponential family distribution assumption for the response variable (y) is relaxed and replaced by a general distribution family, including highly skew and/or kurtotic continuous and discrete distributions. The systematic part of the model is expanded to allow modeling not only the mean (or location) but other parameters of the distribution of y as linear and/or nonlinear parametric and/or additive non-parametric functions of explanatory variables and/or random effects [40]. For our GAMLSS model, it is assumed that independent observations y i , for i = 1, . . ., n have distribution function F Y (y i |θ i ) with θ i = (θ i 1 , . . ., θ i p ) a vector of p distribution parameters accounting for position, scale, and shape. Usually p is less than or equal to four, since lower values provide enough flexibility for most applications. Given an n length vector of the response variable y T = (y 1 , . . ., y n ), let g k (.), for k = 1, . . ., p be monotonic link functions relating the distribution parameters to explanatory variables and random effects through an additive model given by where θ k and η k are vectors of length n, e.g., θ k = (θ 1 k , . . ., θ n k ). β k = β 1 k , . . ., β jk is a parameter vector of length J k and J k itself is the number of explanatory variables for θ i k , X k is a known design matrix of order n×J k , Z jk is a fixed known n×q jk design matrix and γ jk is a q jk -dimensional random variable with γ jk ∼ N q ij (0, G −1 jk ) where G −1 jk is the generalized inverse of a q jk × q jk symmetrix matrix G jk = G jk (λ jk ), which depend on a vector of hyperparameters λ jk . In (2), the linear predictors η k , for k = 1, . . ., p are comprised of a parametric component X k β k (linear functions of explanatory variables), and additive components Z jk γ jk (linear functions of stochastic variables, also denoted as random effects).
The parametric vectors β k and the random effects parameters γ jk , for j = 1, 2, . . ., J k and k = 1, 2, 3, 4 are estimated within the GAMLSS framework (for fixed values of the smoothing hyperparameters λ jk ) by maximizing a penalized likelihood function p given by is the log likelihood function and F is the selected population PDF, and γ jk and γ jk are two independent normal distributions. In this article, we have tested different F by maximizing the penalized likelihood function to find the best fit for the available data. The results of this study are presented in Section IV-B. The selected GAMLSS distribution is the modified sinh-arcsinh (SHASHo) with the following definition: and Here, (σ, μ, ν, τ ) are η k in (2) for k = 1, . . ., 4, which will be estimated based on some explanatory variables [like (7)-(10) in Section IV-B].
The following procedure is used to generate the prices and ATD scenarios.
1) Using one year of historical data for the prices of dayahead energy and FCR-N markets, we select specific number of scenarios for the second-stage of our model. 2) Having the selected second-stage scenarios, we generate the required number of third-stage scenarios using the modified LSTM-based model in Section III-A. 3) Using historical data, we estimate the PDF of ATD for up and down regulation offers by testing different available distributions in GAMLSS and select the one which has the lower AIC and SBC values.

4)
Knowing the day-ahead energy and balancing market prices obtained from step 2) and 3), we generate appropriate ATD for each third-stage scenario using the GAMLSS model extracted from step 3). The results of this section on the proper GAMLSS model selection are elaborated in Section IV-B. From the adaptability point of view, we should state that proposed GAMLSS-based forecasting method for ATDs, it is fully adaptable with any inputs. It takes different inputs such as power generation, electricity market prices, or even consumption level in each area and returns corresponding coefficients, which shows the dependency of ATD over that input [general form in (2)].

A. Case Study Description
To evaluate the effectiveness of the proposed three-stage stochastic optimization, we implement it on a case study of a river in the Northern part of Sweden located in the SE2 bidding zone. It consists of the cascaded HPPs in the Skellefte river. There are 15 cascaded HPPs with a total generation capacity of 1011 MW. The details of each power plant have been presented in Fig. 3. Price information on the day-ahead and balancing energy markets have been extracted from Nordpool market data [39]. The reservoir levels are considered to be 50% of the maximum capacity and should not be depleted below 90% of this value at the end of the operation date. Note that the operator is obliged to set aside some part of its volume in the balancing capacity market with specific minimum requirement announced from TSO. He needs to decide how to split the available capacity in the mFRR capacity, FCR-N, day-ahead market, and mFRR energy markets. We here assume that one actor operates the whole river and this actor then provides one offer-function per hour for the whole river. We also assume k dir = 0.05 and k m = 0.98. All simulations have been carried out in a standard PC with Intel Core i7 CPU and clock rate of 2.1 GHz using no more than 16 GB of RAM. The code is compiled using CPLEX solver in GAMS.

B. Proposed GAMLSS-Based Scenario Generation Algorithm
Based on the algorithm proposed in Section III, the GAMLSS package in R has been employed to find the best distribution fit for the training data. It will linearly regress the parameters of the assumed distribution on the explanatory variables. We have considered the price premium of the regulation markets and day-ahead energy markets (λ pre,t ), wind generation (P wind,t ), and total power production (P pro,t ) in the selected bidding zone with the coefficients h μ,pre , h μ,wind , and h μ,pro , respectively. Following equations show the connection of the explanatory variables to the specifications of the selected distributions. In these equations, μ t is the mean, σ t is the standard deviation, ν t is the skewness level, and τ t is the kurtosis factor In the abovementioned equations, D stands for the direction of the ATD, which is upward or downward. Table II shows the qualification metrics of the different distributions for the training data. The Akaike information criterion and Schwarz Bayesian information criterion (SBC) are two criteria that are asymptotically justified as predictors for the degree of fit in a training data [23]. The lower the values of these two factors, the better the quality of the fitness for the specific distribution. It is clear from the results in Table II that the SHASHo is the best candidate for both upward and downward ATD models as it has  the lower AIC and SBC values. Another reason for this selection is that SHASHo distribution is best suited for the datasets in which there are many zeros or near to zero values (as it is in our case) [40]. This fact is also depicted in Fig. 4, which shows the concentration of the generated data near zero, due to the high level of Kurtosis factor, compared to the normal distribution. Fig. 5 shows the residual of the fitted GAMLSS model to the historical data. The check of the normalized quantile residuals would provide a proper way for testing the adequacy of the fit [40]. Fig. 5 shows the plots of residual: (top left) against the fitted value of mean; (top right) against an index (i.e., number of hours in data); (bottom left) a nonparametric kernel density estimate; (bottom right) a normal Q-Q plot. It is expected that the fitted residuals behave approximately as normally distributed variables (even though the original observations are not necessarily normal); thus, a normal Q-Q plot of the residuals is appropriate. The statistical moments of the quantile residuals are stated in Table III, which shows proper fit as mean and variance are close to zero and one, respectively.
Tables IV and V show the intercepts and coefficients obtained from GAMLSS for the different explanatory variables in the upward and downward trends. These tables show that the ATD has the most dependencies on wind power generation data compared   to the price premium and total production. It can be explained by this fact that with increasing the penetration level of wind power in the grid, its uncertain electricity generation changes the TSO's decisions on the activation level of balancing energy offers; hence, the ATD will be affected. Fig. 6 demonstrates the time-varying distribution percentiles of the ATD for up and down regulations in 100 scenarios shown as shaded bands around a central mean line.

C. Revenue Breakdown Analyses of the Proposed ATD
In the previous studies (e.g., [15], [17], [20]), which have not considered ATD of balancing energy offers, there is an illegal arbitrage between day-ahead and balancing energy markets. On those works (considering high prediction for the value of stored water), one possible strategy is to allocate some (or full) part of available capacity to the day-ahead energy market and then buy it back in the downward balancing energy market with a lower price without any electricity generation obligation. This act is a clear sign of increase-decrease game, which has been observed in many parts of the world [21]. To alleviate this problem, we have considered the following three strategies for the ATD. 1) Naive case: In this case, we have full ATD, which means that all the submitted offers in the mFRR energy markets will be activated for one hour. This assumption provides the possibility for the operator to offer with full capacity to the up regulation market instead of day-ahead or FCR-N market, which provides higher revenue while it is not realistic. There reason is that not all of up regulation balancing energy offers would be activated and based on the imbalance level, some part of it will be activated by TSO.
2) Proportional ATD: In the proportional ATD, we have used the known data to find an estimate of ATD. These data have been calculated using the total upward (downward) activated volume of all offers in a specific biding zone, in which our case study is located, divided by the total procured upward (downward) volume in that hour. We solve this optimization problem before the day of operation while the information we have used in the proportional ATD calculation will be revealed after the day of the operation; however, it is closer to the reality than the naive case. This calculation is useful to know what would be the expected revenue breakdown if we had known ATDs.
3) Proposed GAMLSS-Based ATD: In this strategy, we have estimated the ATD based on the approach proposed in Section III-B. The results are stated in Tables VI and VII as they show two cases of expected value of stored water in (1a). As it is clear from the results, the revenue stream in the proposed ATD has shifted toward the day-ahead market compared to the naive case because it is not possible to be remunerated with full ATD in the balancing markets. In fact, TSO will activate the accepted offers for specific time duration. Although the total revenue is lower than the other cases, but it is more realistic to expect this revenue as we take into account the possibility of not being dispatched.
To clarify the previous case studies, we can state that we have extracted one year of historical data to generate scenarios for the balancing market prices in the second stage of our problem formulation and LSTM-based scenarios generation technique in the third-stage for mFRR energy market prices. Here, we have assumed uniform distribution for the scenarios in both cases. But for the ATDs in both directions, we have two approaches based on the proportional ATD and GAMLSS-based ATD as we have described in this section. In the proportional ATD approach, we have used the uniform distribution of the scenarios calculated based on the description presented in the related paragraph in this section. But, for the GAMLSS-based ATD, we have used GAMLSS to find the best distribution fitted to the input data [which we have used those historical data to generate statistical moments of the selected distribution, i.e., SHASHo in (7)-(10)] and then used this distribution to generate scenarios.

D. Profitability Assessment of FCR-N Market Participation
In this section, the profitability of participating in the FCR-N market is explained. Table VIII shows the revenue breakdown from different markets with and without considering participation in FCR-N market and Table IX shows their allocated volume. As it shows, the total revenue is increased as more volume has been allocated in day-ahead energy and FCR-N markets. Participation in the FCR-N market is more beneficial than the mFRR capacity market due to the higher prices, especially at night, and no obligations in the real-time market. The optimal offer-functions are calculated and presented in Fig. 7. This figure shows the optimal offer-functions submitted by HPP operator for participation in the day-ahead energy market. To illustrate the trend, the following two scenarios have been considered. 1) High FCR-N Price Hours: During night hours, the FCR-N prices are high due to the lower available capacity in the network for down regulation and the need for bidirectional flexibility from FCR-N providers. Thus, HPP operator wants to maximize its participation in this market to gain more profit. Due to the less obligation compared to the mFRR capacity market, he prefers to offer more in FCR-N market and less in the mFRR capacity markets, and considering (1h), it is better to be more dispatched. Therefore, the offer volume increases as it can be seen in Fig. 7 at t = 4 and t = 7. 2) Low FCR-N Price Hours: During the day, the FCR prices decrease due to the more availability of capacity to be dispatched through market participants for the frequency fluctuations in the grid. Thus, the HPP operator seeks to allocate more volume into the mFRR capacity markets and less to the FCR-N market. As a consequence, the mFRR energy market in the upward regulation becomes more profitable than the day-ahead energy market and decreases the offer volume of day-ahead energy market as shown in Fig. 7 at t = 15 and t = 21.

E. Evaluation Metrics of Proposed Scenario Generation Methods
Energy score (ES) is a skill score that evaluates the quality of generated scenario sets [41]. It provides a direct comparison between scenario sets. This is a negatively oriented value, i.e., it is inversely related to the skill of the scenario set. Equation (11) shows how we can calculate this score: ξ ω − ξ ω (11) where . denotes the Euclidean norm, y is the realization of the random process under study, ξ ω and ξ ω are independent scenarios generated by the same forecast technique, and N Ω is the total number of scenarios. In Fig. 8, we compare the results of  the forecasted up regulation prices for 120 h with the available historical data and their corresponding error. It demonstrates proper fit with total root mean square error (RMSE) equal to the 53.1 € /MWh. Fig. 9 shows the ES for two sets of data for the mFRR energy market in upward and downward trends. It compares the EC score for the following three approaches: conventional LSTM (C-LSTM), modified LSTM (M-LSTM), and SARIMA. The difference between C-LSTM and M-LSTM is that in M-LSTM, we have the updated forecasting method based on the available information as proposed in Section III-A. With a lower ES at each instance, LSTM methods can perform better and generate scenario sets with higher skill than other statistical approach examined in this article. The relative out-performance of LSTM methods can be attributed to their internal memory capacity and the ability to store information over multiple time steps.

V. CONCLUSION
This article proposes a multistage stochastic program to find the optimal multimarket offer-function for the participation of cascaded HPPs in different electricity market platforms. The day-ahead, FCR-N, and mFRR markets are considered in our article. In the multimarket offer-function calculation, the ATD is a determining uncertain factor. We propose a modified GAMLSSbased machine learning approach for modeling the uncertain ATD in the offer-function calculation. The ATD consideration is an important factor toward a more realistic multimarket offerfunction calculation for which full-hour activation of balancing offers is a naive assumption. Furthermore, we show that introducing a modified LSTM approach can improve the scenario generation quality for balancing energy market prices. Also, the inclusion of FCR-N market with the mFRR capacity market changes the offer-function of the HPP operator in the day-ahead energy market. Future work can use machine-learning techniques to model the functional relationship between explanatory and response variables in our proposed GAMLSS-based distributional regression model. Also, our developed statistical tool can also be used for estimating and predicting the economic welfare in the electricity markets. The inclusion of risk assessment techniques, such as conditional value at risk (CVaR), would be a valuable study to investigate the revenue at risk when comparing different strategies for ATDs.