Intraday renewable electricity trading: advanced modeling and numerical optimal control

As an extension of (Progress in industrial mathematics at ECMI 2018, pp. 469–475, 2019), this paper is concerned with a new mathematical model for intraday electricity trading involving both renewable and conventional generation. The model allows to incorporate market data e.g. for half-spread and immediate price impact. The optimal trading and generation strategy of an agent is derived as the viscosity solution of a second-order Hamilton–Jacobi–Bellman (HJB) equation for which no closed-form solution can be given. We construct a numerical approximation allowing us to use continuous input data. Numerical results for a portfolio consisting of three conventional units and wind power are provided.


Introduction
To counter global climate change, renewable power sources substituted fossil fuel plants and provide now a substantial part of the electricity production. Due to the intermittency of renewable power, short-term electricity contracts have gained importance on electricity exchanges such as the European Power Exchange (EPEX Spot a ). In particular, continuous intraday trading, which allows trading of contracts until 30 minutes before delivery, is used to respond to short-term changes. The trading volume within the German intraday market (IDM) area increased from 26 TWh in 2014 to more than 50 TWh in 2018.
A similar trend has been observed in other market areas or countries in which some markets or sub segments, e.g., continuous trading of hourly products in Belgium (since July 2018) or the 30-min continuous trading in Germany and France (since April 2017), have been developed. Another instrument for integrating renewable energy markets is the Xbid project, which aims at establishing a common pan-European continuous intraday market to strengthen liquidity. All these developments trigger a need for mathematical modeling of such trading as a basis for deeper understanding, optimization and control. Moreover, a mathematical model is the basis for numerical simulations.
However, in contrast to day-ahead market (DAM) modeling, literature dedicated to the continuous intraday electricity markets is scarce. In particular, appropriate mathematical models describing the main characteristics of these short-term markets have not been developed so far. Most important price drivers have been identified in [10] explaining about 75% of the price variance, [20]. Moreover, strong statistical evidence was found that information on the specific fundamental factors significantly affect the intraday prices, which appear as transaction prices within a trading period, [13,14].
Concerning mathematical modeling early work in the field of integration of renewables into short-term markets using stochastic optimization can be found in [18]. The minimization of incurred intraday costs of wind producers while maintaining the balance of production forecast and sales is considered in [11]. This short-term trading model also considers the impact of the wind producer on prices without intraday price uncertainty, however. A discrete time decision model with intraday prices following a geometric Brownian model and wind production error forecast following an arithmetic Brownian motion has been introduced in [7,8]. In that framework, the power producer is supposed to have no impact on intraday prices.
To overcome the aforementioned weaknesses, we have been inspired by [1] to model the continuous intraday market for electricity. In [1], a Hamilton-Jacobi-Bellman (HJB) equation was derived for determining an optimal trading strategy by modeling the dynamics of the electricity market by stochastic differential equations (SDEs) and formulating a corresponding value function to be optimized. The specific market model allows one to solve the arising HJB analytically, i.e., the authors derive a solution formula.
The starting point of this paper (which builds upon and extends [9]) is a statistical analysis of EPEX SPOT data, which shows that some of the model assumptions in [1] are in fact not satisfied under real market conditions. Thus, we introduce a more sophisticated model. The arising HJB equation can no longer be solved analytically; the value function is shown to be the unique viscosity solution of this HJB equation. Thus, we need an appropriate numerical scheme. From an economical point of view, the main new ingredients of our model are: 1. Portfolio of renewable and conventional energy represented by a cost function that reflects the stepwise merit order of a portfolio rather than a systemwide quadratic function; 2. Pricing model using time-varying half-spread and being capable of representing time-varying liquidity; 3. Approximation of market data for half-spread and instantaneous price impact; 4. Variable penalty depending on the state of the market at final time; 5. The model relies on data which are observable on the market (see Sect. 3). The main focus of this paper is a novel application-related modeling of the intraday trading and the determination of a numerical approximation for this problem. We show an example of a real-world problem, compute the optimal trading strategy and investigate the impact of various involved parameters. The remainder of this paper is as follows: In Sect. 2, we introduce the new model and the arising HJB equation and Sect. 3 details the involved parameters, which we obtained from empirical analyses. In Sect. 4, we describe the numerical method for determining the viscosity solution of the arising HJB. More-over, we present corresponding numerical experiments. We finish by conclusions and an outlook in Sect. 5.

A new mathematical model
In order to take both renewable and conventional generation into account, our model is based upon the consideration of an agent owning both kinds of power plants and aiming at selling a combination of renewable b and conventionally produced electricity. In detail, depending on the weather forecast and the expected price at the final time, such a combination of electricity is sold at the day-ahead market (DAM). With the result of this initial trading, the agent starts to act on the continuous intraday trading market (IDM) aiming at maximizing her profit by determining an optimal trading strategy as well as an optimal production strategy of conventional power. c All quantities entering the model are described below.

Day-ahead and intraday electricity trading
Consider a delivery hour h on day d. The day before, d -1, the day-ahead auction takes place with gate closure at 12 noon. In this auction, each participant can offer (ask) or request (bid) a certain demand of electricity at a specific price. Then, a clearing price is set and power is exchanged accordingly. Next, the continuous intraday trading starts at 3 p.m. on day d -1 and closes half an hour before the actual delivery hour h, see

Dynamics of the electricity market
The dynamics of the market includes the forecasted renewable power production and the price process. The latter one is influenced by the current trading activity of the agent. We use stochastic processes and derive stochastic differential equations (SDEs), see [19] for background.

Renewable production forecast
By D = (D t ) 0≤t≤T we denote the forecasted production of renewable produced electricity during the trading session. Following the idea of [15,22], we assume that forecast updates are a consequence of new information and hence lead to random changes in t. Therefore, the uncertainty is modeled by means of the dynamics where σ D is the volatility and (W t,D ) 0≤t≤T is a standard Brownian motion. For the sake of simplicity, this variable is unbounded, whereas in the real world, there are restrictions by zero (no wind) and the maximum capacity of the wind farm. We denote by D t,d the solution of the SDE starting from d (the current capacity) at t.

Agent's position
The financial position resulting from the agent's trading activity is denoted by X = (X t ) 0≤t≤T . The agent participates in the intraday market (IDM) with continuous trading at rate q t ∈ Q := (q min , q max ) ⊂ R, where q t > 0 implies actively buying and q t < 0 implies actively selling. d The dynamics of her financial position are then and we denote by X q;t,x the solution of the SDE starting from x at t depending on the current trading rate q t . For t = 0, x 0 is the amount of electricity sold on the day-ahead market.

Permanently impacted mid price
First, we denote by Y = (Y t ) 0≤t≤T the permanently impacted mid price, i.e., the sum of the mid price of energy and the permanent impact of the agent's trading, the latter one modeled by some function ψ : R → R. The dynamics of Y is modeled by the SDE where μ Y is the drift, σ Y is the volatility and (W t,Y ) 0≤t≤T is a standard Brownian motion. We denote by Y q;t,y the solution of the SDE starting from y at time t depending on q t .

Transaction costs
When actively buying (selling) on the IDM, there are costs in addition to the permanent impact. Those are referred to as transaction costs and will be incorporated in our model.

Execution price
The execution price is the price the agent pays (receives) when actively buying (selling). We require a more advanced approach of the pricing model as in [1], where the half-spread (see below) and its time variability as well as the time variability of the execution costs are ignored. Incorporating these effects, the execution price depends on several quantities to be introduced now. The half-spread is defined as the half of difference of the best ask and the best bid price. This data is observable on the market. While in reality the bid-ask spread and hence also the half-spread is stochastic, we model it in terms of a deterministic function h : [0, T] → R + , which reflects the typical shape over the trading period.
The other component of our model for the execution price is the execution cost reflecting the costs that are incurred due do executing limit orders with prices worse than the best price when buying (selling) actively. While they are also stochastic in reality, they will be described in terms of a deterministic function ϕ : [0, T] × (q min , q max ) → R which reflects the typical shape over the trading period.
With all these quantities at hand, the execution price P q;t,y = (P q;t,y s ) 0≤s≤T is modeled as where (as usual) Hence, (4) means that the execution price is the permanently impacted mid price plus (minus) the time-varying half-spread and plus the time-varying execution costs, which extends the model in [1].

The profit: the objective function for optimization
It is the aim of the agent to maximize her profit consisting of the running and the terminal profit explained next. This will later serve as the objective function for the arising optimization problem.

Running profit
The total running profit from the continuous trading in the intraday market is given by with the trading rate q = (q t ) 0≤t≤T .

Terminal profit
The profit gained at the end of the trading period consists of scheduling conventional power production which incurs costs and placing a final market order.
Conventionally produced energy At the end of the trading session, the agent chooses how much electricity ξ ∈ R + 0 she will produce during the delivery period. For doing so, she has n conventional units available with each being able to operate between their specific minimum κ min i and maximum generation capacity κ max i , i = 1, . . . , n. The marginal costs c i of unit i, i.e., the costs for producing 1 MWh of electricity, are assumed to be time-independent (constant). The chosen strategy of the agent thus consists of deciding to activate or deactivate unit i (modeled by a binary variable a i ∈ {0, 1}) and choosing the respective amount ] for each unit. Thus, the resulting amount of produced electricity is A straightforward strategy would be that the agent will activate her units in ascending order of marginal costs, starting with the 'cheapest' one. The arising total cost of power production then reads which is a piecewise linear but discontinuous function, the derivative of which is a sum of piecewise constants and a number of Delta distributions.
Final market order Furthermore, the agent also has the option to place a final buy or sell market order. In fact, obtaining values for D T , X T and Y T , the agent optimally tries to reach her desired demand for power by buying/selling the amount ξ + X T + D T ∈ R, i.e., what has been traded already minus conventional and renewable production, by a final market order. The costs associated with the final market order are the costs due to crossing the half-spread h(T) and the costs due to potentially executing limit orders with prices worse than the best bid/ask price. The latter costs are modeled by the function α : R → R.
Terminal payoff It turns out to be convenient to introduce the variable Z t := X t + D t , i.e., the sum of the forecasted production from renewables and what has been sold by the agent so far. This, in particular reduces the dimension of the optimization problem to be solved numerically, which is crucial for efficiency. The terminal payoff/profit then reads

The resulting optimization problem
Now, we have all ingredients at hand to formulate the optimization problem in terms of a HJB equation.

Value function
The value function corresponds to the agent's cash, so that an optimal strategy yields maximal cash. Denoting by Z q;t,z the solution of the SDE dZ t = dD t + dX t = q t dt + σ D dW t,D starting from z at t, the resulting value function V : [0, T] × U → R reads for Q := (q min , q max ), see, e.g. [24,Ch. 4], where U := Y × Z ⊂ R 2 is a closed convex set (usually a rectangle, in order to ensure wellposedness of the optimization in (9), [6]). In fact, from our above modeling, we would obtain Y = Z = R (range of values for the execution price and energy production). However, without additional conditions, the optimization problem is not well-posed on such unbounded domains. Moreover, the numerical solution using standard discretization techniques is-at least-not straightforward on unbounded domains and would require sophisticated schemes, see e.g. [12]. In order to overcome these difficulties, we define U in terms of closed sets (intervals) Y and Z and cut-off the problem on their boundaries. This requires to set boundary conditions to ensure well-posedness to be explained below.

The Hamilton-Jacobi-Bellman (HJB) equation
Following the well-known dynamic programming principle (e.g. [24,Ch. 4]), we derive the HJB equation: for (t, y, z) ∈ [0, T) × U with terminal condition V (T, y, z) = g(T, y, z) for all (y, z) ∈ U . Note, that in (10) we used the notation q(t) instead of the previously used q t for the following reason: In the numerical realization, we treat q as a function q(t) and not as stochastic process q t . One can show that (10) with terminal and Dirichlet boundary conditions is well-posed and that the unique viscosity solution is the value function in (9), [5]. Due to the form of (10), we cannot expect a first-order condition for the control q and we have to resort to numerical solvers.

Boundary conditions
As already mentioned above, we need to truncate the HJB to a bounded domain U , which requires to prescribe appropriate boundary conditions on ∂U . Such an approach is wellknown also from numerical option pricing when solving the Black-Scholes (BS) equation.
In that case, the definition of corresponding Dirichlet boundary conditions is canonical since small errors in the boundary values cause only small errors in the solution on all of U (i.e., the BS equation is stable w.r.t. boundary values). The situation for the HJB equation is fundamentally different as small changes in the Dirichlet data immediately severely change the solution on the whole domain. Hence, the definition of appropriate boundary conditions is a delicate task.
We are going to describe our approach to prescribe appropriate boundary conditions. Our point of departure is the before-mentioned IDM-model in [1], which uses a much simpler HJB as the one derived above. This causes the fact that the HJB in [1] allows for a closed formula for the solution. We consider such a somehow simpler HJB, solve it explicitly and use the boundary values of that HJB as Dirichlet data for our more sophisticated HJB. Of course, one could also use other simplified models as long as the resulting boundary conditions turn out to be meaningful. To be specific, we consider a simpler HJB model consisting of the following ingredients and simplifications: 1. By omitting the sign-part in (10) and replacing it by the positive sign, the function within the supremum is differentiable. 2. We replace h(t) byh := h(T), i.e., the half-spread at the terminal time. 3. The next simplification concerns the temporary price impact ϕ. The most simple situation would be a linear approximation, e.g., ϕ(t, q) := kq, i.e., stationary and linear in q. Here, k ∈ R + is a constant, which has been derived by approximating our mid price data with a second order polynomial p and then setting k := p(T). 4. The permanent price impact vanishes, i.e., ψ ≡ 0. Thus, the simplified HJB now takes the form In (11), we can explicitly determine the supremum by finding the root of the first-order derivative w.r.t. q of the term in {· · · }. The corresponding first-order necessary conditions yield the optimal control as follows (recall from 3. that k > 0) Inserting this optimal control into (11) yields the following PDE for the unknown v as a function of (t, y, z) To solve (12), we make the following polynomial ansatz i.e., a polynomial of degree 2 with the coefficient functions a i : [0, T] → R to be determined. Plugging this form of v into (12) yields a system of six Riccati equations for the unknown functions a i , i.e., In order to solve this system of ordinary differential equations, we need initial conditions for the functions a i . Due to the arguments (Tt) in all those functions, the desired initial conditions boil down to terminal conditions of v, i.e., the function g in (8). Recall, that g in particular contains the function C in (7), which is piecewise polynomial but discontinuous (see, e.g. Fig. 6 for an example of such a function-clearly exhibiting jumps), which prohibits a closed form solution of the Riccati system. Thus, we use a least squares approximation of g in terms of the above polynomial v(T, y, z). Doing so, we obtain initial values for the above mentioned functions, say a i (0) = a i,0 . With these values at hand, we solve the initial-value problem of the Riccati system by Maple™ and obtain v. The boundary values of v w.r.t. the variables y and z are then used as Dirichlet conditions for (10).

Data
Our model described above relies on several parameters, which we summarize in Table 1.
In this section, we describe how this data can be obtained from market observations and empirical data analysis.

Generation portfolio
As mentioned before, the agent's portfolio consists of renewable and conventional generation capacity. We describe how to retrieve realistic market data.

Renewable electricity standard (RES) portfolio
We equip the sample agent's portfolio with 500 MW of aggregated renewable generation capacity. As mentioned before, for the sake of simplicity and due to limited availability of historical data, we restrict the aggregated capacity to arise solely from wind farms. Furthermore, following [17], we assume that the wind farms in the portfolio are dissimilarly located within the considered hypothetical market area. This last assumption allows us to estimate the parameter σ D for D t in (1) using aggregated forecast data, which is, in contrast to site-specific data, publicly available. Specifically, we choose hourly wind power forecasts for the French market area provided by [21]. To be able to transfer the characteristics of the historical data set to our current application, we first normalize all forecasts on the average installed capacity per month. We then determine all updates between two adjacent forecasts of the same forecast path. Finally, we observe a volatility of approximately 0.01 per installed MW and hour (2016: 0.008, 2017: 0.008, 2018: 0.010) in the data. With respect to the renewable generation capacity of 500 MW, we therefore set σ D = 5 MW. e Finally, we choose D 0 also from publicly available data.

Conventional generation
For the agent's conventional portfolio, we consider n = 3 units, namely a hard coal fired plant, a combined cycle gas turbine (CCGT) and an open cycle gas turbine (OCGT), with the parameters shown in Table 2. The marginal costs c i of each unit represent idealized values for the respective technology class. They also consider that an increase in flexibility-here the reduction of the so-called deadband between zero and production at minimal capacity κ min i -reduces the efficiency of the unit. Moreover, we assume that the start-up decision a i ∈ {0, 1} of a unit does not require a lead time.

Mid price drift
We assume that the drift of the mid price on the IDM consists of two parts as implied by (3), namely mid price changes due to time evolution on one hand and mid price changes due to agent's trading (causing irreversible price impacts) on the other. Building upon [4], we perform an empirical analysis of these two parts. In particular, we study price changes over longer periods of time as well as their relation to net order flow (i.e., buy minus sell). While in that study mid price changes and net order flow are considered over periods of five minutes, we consider the differences between the day-ahead market (DAM) prices and volume-weighted average IDM prices f as well as the net order flow over the entire trading window. While we conjecture the DAM prices to be close to the mid prices after market opening, the volume-weighted average IDM prices are usually different from the mid price before end of trading. Nevertheless, we prefer mid prices as we reckon that they better reveal the evolution of the price over a longer period of time and also mirror the relation of this evolution to net order flow. Of course, the observed price changes appear between the beginning and the end of the IDM trading. As a consequence, the data does not show whether either of the components typically changes over the trading window. Concerning the (deterministic) dependence of the price on time, we simply choose a linear dependence, i.e., constant drift as μ Y (t) ≡ μ Y .
Next, we assume a linear relation of order flow and price change, which is at least not contradicted by the scatter plot in Fig. 2. Thus, we assume ψ(q t ) = bq t for the permanent price impact with a constant b ∈ R.
We obtain estimates for b and μ Y from least squares fits to the data. For b, we get 0.0017e/MWh 2 and significance at a 0.1% level, indicating that the net order flow has a positive impact on the price change. The drift μ Y is obtained as 0.0433e/MWh per hour and to be significant at a 1% level, indicating that the price slightly tends to increase over time. Concerning the sign of the drift, we find mixed evidence in the literature, [10,14].

Transaction costs
As we pointed out earlier, it is a major difference between our model and previous research that we also include transaction costs. The data analysis on which their modeling is based is presented below. We recall that transaction costs include execution costs modeled by the function ϕ and the half-spread h, see (4).

Data
We use order book data from the EPEX SPOT-operated market for hourly delivery contracts with Germany/Austria g as delivery area from the second quarter of 2016 (referred to as Q2/2016 in the following) to empirically analyze the transaction costs mentioned in Sect. 2.2. The dataset comprises (i) all orders with Germany or Austria as delivery area which entered into the order book, (ii) all orders with Germany or Austria as delivery area which caused execution of an order resting in the order book with delivery area other than Germany or Austria and (iii) all orders with delivery area other than Germany or Austria which caused execution of an order in the order book with Germany or Austria as delivery area.
Hence, not all orders in the order book for the German/Austrian delivery area which were visible for market participants are contained in the dataset. Based upon these data, we are now going to describe how we obtained values for the quantities entering the transaction costs.

Half-spread
In identifying typical half-spread functions h entering the optimization problem, we build on research on bid-ask spreads (BAS) on the NYSE h stock market. In [16], the pattern of BAS of NYSE stocks over a trading day is analyzed. To this end, the authors divide each trading day in their sample into one-minute intervals and compute for each interval and stock what they refer to as the time-weighted BAS to be explained next. Consider some time interval I i := (T i-1 , T i ] ⊂ [0, T] and assume that the BAS changes N i times at ) by BAS j , j = 0, . . . , N i setting t (i) 0 := T i-1 and t (i) N i +1 := T i . Then, the time-weighted BAS on that interval I i , denoted by BAS i , is defined as Then, [16] suggests to determine BAS i for each interval of the trading day and for all stocks in their sample. In this paper, we adopt the approach in [16] and divide the trading period into intervals of 5-minute length. Furthermore, we only consider the hourly delivery contracts with delivery starting at 12 noon in Q2/2016. The resulting data is visualized in Fig. 3 as a two-dimensional histogram of all the BAS i in the sample for the last 17.5 hours of trading. The blue line reflects the means of the BAS i in each interval. We observe a significant decrease of the mean time-weighted BAS from ≈8e to ≈5e/MWh at the beginning of the time window and a subsequent nearly constant behavior for about ten hours. Given a tick size of 0.1e/MWh (Q2/2016), it is remarkable that this plateau is fifty times higher than the tick size. Five hours before the end of the trading, the mean time-weighted BAS decreases to ≈1e/MWh followed by a sharp increase to reach the final ≈2e/MWh. The pattern over the last five trading hours is quite similar to the crude reverse J shape reported in [16] for NYSE stocks. In order to model the temporal behavior of the BAS i , we employ a polynomial of degree seven according to the Akaike information criterion (AIC) [2], which is depicted by the red line in Fig. 3. This polynomial is the half-spread function h mentioned in the previous section. Clearly, this has a significant smoothing effect and could e.g. be replaced by other approximations as well.

Execution costs
The execution price in (4) depends on the typical half-spread and the typical execution costs, which reflect the negative impact on the price realized by a market participant when buying or selling with market orders. For the empirical analysis of the execution costs we need to consider the order book. Similar to the analysis for the half-spread, we start by determining time-weighted prices and volumes over 5-minute intervals on the different order book levels. i Given some time-weighted price on an order book level, it may occur that the time-weighted price on a lower level is better. Therefore, we sort prices and volumes in descending/ascending order on the buy/sell side. This approach requires the availability of the entire order book over the trading period. Otherwise, missing data techniques could possibly be used.
We assume a linear relationship between trading rate and execution costs, with the parameter (function) k(t). For estimating k(t), we build upon [4] and analyze how the order books absorb market orders of different sizes. To this end, we consider market orders with volume 1, . . . , 200 MWh for each interval and market side. Then, we collect those order book levels required to cover the volume of the market order. We multiply the order book level prices in that collection by their volumes, sum them up and divide by the volume of the market order. From the resulting price we subtract the best price on the same market side to obtain the respective market order's execution costs. Then, we fit a linear model by least squares to obtain k(t). Considering again hourly delivery contracts with delivery start at 12 noon in Q2/2016, Fig. 4 shows the obtained values for k(t) in the form of two-dimensional histograms as well as means per time interval (blue line). Similar to the average half-spread, the average execution costs exhibit a decline after market opening to 0.025-0.05e/MWh 2 . After a rather stable ≈10 hours period, they further decline to ≈0.01e/MWh 2 . We observe a slight increase just before the end of trading.
Hence, the shape of the average execution costs is rather similar to that of the average halfspread. We model the typical temporal behavior of the execution costs with a polynomial of degree six for the buy side and degree eleven for the sell side (according to the AIC).
Remark 1 Note that the above approach is not compatible with the model in (13). For determining a model for the execution costs, we compared market orders and the order book. However, the agent's action in our model is trading at some rate. Of course, in reality, market participants will not act by trading at some rate, but merely actively place market orders by some strategy. This means, we cannot observe how different rates enter the order book. On the other hand, the above approach is mainly used to calibrate parameters for our model.
A possible strategy to remedy this shortcoming could be to determine a relation between execution costs in some time interval and a constant trading rate (instead of the volume itself ). For example, for execution costs evolving from a 1 MWh market order, a trading rate of 1/5 MWh per minute would be required to yield that volume after 5 minutes of trading.

Terminal order book
At the end of the trading period, we assume that the agent liquidates remaining inventory by means of a final market order (instead of letting the volume run into the balancing market as is done in [1]). This means that we need a typical order book at the end of trading for calibrating our model.
To this end, we consider the buy and sell side separately and determine for each contract the difference between the prices on the different order book levels and the best price. Then, we average all these price differences and volumes on the same order book level. The results are shown in Fig. 5.
Given that the average volumes associated with the best bid and ask are around 16 MWh, typically there is still some volume at the end of trading which can be sold/bought at zero execution costs. The volumes associated with the order book levels beyond the best-price level slightly increase. While the average price differences associated with the best bid and ask are obviously 0e/MWh, in absolute terms they are ≈1e/MWh for the first level, ≈5e/MWh on the fifth level and 20-30e/MWh on the tenth level. Hence, prices beyond the best level obviously worsen quite substantially.
We consider both the typical order book on the buy and sell side at the end of trading to specify the penalty function α. Recall, that δξ := ξ + X T + D T denotes the 'untraded' amount (for which a penalty needs to be paid) with ξ being defined by (6). Furthermore, we only consider the first L levels of both the buy and sell order book and truncate the volume on the last level (L) such that the overall volume is 100 MWh on both market sides. Let δp buy ≤ 0 (δp sell ≥ 0) denote the difference between the price on the -th buy (sell) order book level and the price on level zero, = 1, . . . , L. Furthermore, let λ buy , λ sell be the maximum volume available on respective side of the -th order book level. Then, α is defined as It is easily seen that α is continuous at δξ = 0. The sign in (14) results from (8) and the motivation that the penalty α should lower the agent's profit. This implies that α should be positive if the agent needs to buy and negative if she needs to sell at the end of the trading period.

Numerical solution of the HJB equation
In this section, we describe our numerical method for (approximately) solving the arising HJB. Moreover, we report on results of a sample numerical experiment concerning (10) using the following data We use a finite difference discretization from [23] with 56 × 301 points in space and 100 points in time. In particular, central differences are used for the approximation of the first-order terms with additional artificial diffusion, which results in a stable, consistent and monotone scheme converging to the viscosity solution, [23]. We use the well-known policy iteration [3] in every time-step and the control is maximized over a discrete set (as no first-order conditions are available). Finally, the optimal conventional generation is computed as the maximum value of (8) w.r.t. ξ using Matlab's j intlinprog with the interior point method.
Our results for the optimal conventional generation ξ are displayed in Fig. 6. Let us comment on the case where Z T = -500 MWh. As long as the final mid price is below 25e/MWh, the agents buys the maximal amount of 145 MWh (recall, that y ∈ [-1645, 145]) and uses the power plant with the lowest marginal costs (hard coal) accordingly, i.e., the remaining 355 MWh. Once the final mid price is 25 to 35e/MWh (i.e., above the marginal cost of hard coal, but below the marginal cost of CCGT) it is optimal to produce at maximum capacity with the cheapest conventional power plant (i.e., 500 MWh by hard coal) and no final market order is required. If the final mid price exceeds 35e/MWh, the agent sells as much electricity as possible (145 MWh) and produces exactly that amount with the CCGT plant at 35e/MWh, which is possible because its capacity is 100-400 MW. Finally, no matter how high the final mid price is, the OCGT unit with the highest marginal cost is not used, since there is not enough sell volume on the market. These results are clearly reasonable. Figure 7 shows the optimal trading rate over the trading window t ∈ [0 h, 17.5 h]. In both cases, we fix Z t ≡ -499.4 MWh (the non-integer numbers arise from the discretization w.r.t. y and z). For the mid price, we choose Y t ≡ 59.25e/MWh (left) and Y t ≡ 13.98e/MWh (right). In the left plot, the trading rate is negative (selling), which is reasonable since Z t ≡ -499.4 MWh means that the agent has only marketed the cheapest power plant and Y t ≡ 59.25e/MWh means that the execution price is above the marginal costs of the second cheapest power plant. Note, that the absolute value of the trading rate substantially increases around 15 h, since half-spread and immediate price impact are minimal there. In the right plot, the execution price is below the marginal costs of the cheapest power plant, the agent buys electricity and reduces the production of the marketed power plant.

Conclusion and outlook
We have introduced an extended model for the intraday market of renewable electricity. As opposed to earlier research, our more sophisticated approach does not allow for a closed solution formula for the desired optimal trading strategy as a function of time. We thus used a numerical scheme for approximately solving the arising Hamilton-Jacobi-Bellman (HJB) equation. The parameters within the HJB equation are market data which we showed to be available by an empirical analysis.
The availability of a numerical approximation scheme allows us now to extend our model to all market participants, so that regulatory constraints can be determined e.g. for reaching desired environmental goals. Moreover, we will use our scheme to further investigate optimal strategies within economically particular relevant market scenarios.