An optimal trading problem in intraday electricity markets

Abstract

We consider the problem of optimal trading for a power producer in the context of intraday electricity markets. The aim is to minimize the imbalance cost induced by the random residual demand in electricity, i.e. the consumption from the clients minus the production from renewable energy. For a simple linear price impact model and a quadratic criterion, we explicitly obtain approximate optimal strategies in the intraday market and thermal power generation, and exhibit some remarkable properties of the trading rate. Furthermore, we study the case when there are jumps on the demand forecast and on the intraday price, typically due to error in the prediction of wind power generation. Finally, we solve the problem when taking into account delay constraints in thermal power production.

Appendix

1.1 Proof of Theorem 3.1

The Hamilton-Jacobi-Bellman (HJB) equation arising from the dynamic programming associated to the stochastic control problem (3.4) is:

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\partial {\tilde{v}}}{\partial t} + \inf _{q \in \mathbb {R}} \left[ q \displaystyle \frac{\partial {\tilde{v}}}{\partial x} + \nu q \displaystyle \frac{\partial {\tilde{v}}}{\partial y} + \mu \displaystyle \frac{\partial {\tilde{v}}}{\partial d} + \frac{1}{2} \sigma _0^2 \displaystyle \frac{\partial ^2 {\tilde{v}}}{\partial y^2} + \frac{1}{2} \sigma _d^2 \displaystyle \frac{\partial ^2 {\tilde{v}}}{\partial d^2} + \rho \sigma _0\sigma _d \displaystyle \frac{\partial ^2 {\tilde{v}}}{\partial y \partial d} + q(y+\gamma q) \right] = 0, \nonumber \\ {\tilde{v}}(T,x,y,d) = {\tilde{C}}(d-x) = \frac{1}{2} r(\eta ,\beta ) (d-x)^2. \end{array}\right. } \end{aligned}$$

The argmin in HJB is attained for

$$\begin{aligned} {\tilde{q}}(t,x,y,d)= & {} - \frac{1}{2\gamma }\left[ \displaystyle \frac{\partial {\tilde{v}}}{\partial x} + \nu \displaystyle \frac{\partial {\tilde{v}}}{\partial y} + y \right] , \end{aligned}$$

and the HJB equation is rewritten as:

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\partial {\tilde{v}}}{\partial t} + \mu \displaystyle \frac{\partial {\tilde{v}}}{\partial d} + \frac{1}{2} \sigma _0^2 \displaystyle \frac{\partial ^2 {\tilde{v}}}{\partial y^2} + \frac{1}{2} \sigma _d^2 \displaystyle \frac{\partial ^2 {\tilde{v}}}{\partial d^2} + \rho \sigma _0\sigma _d \displaystyle \frac{\partial ^2 {\tilde{v}}}{\partial y \partial d} - \frac{1}{4\gamma }\left[ \displaystyle \frac{\partial {\tilde{v}}}{\partial x} + \nu \displaystyle \frac{\partial {\tilde{v}}}{\partial y} + y \right] ^2 = 0, \\ {\tilde{v}}(T,x,y,d) = \frac{1}{2} r(\eta ,\beta ) (d-x)^2. \end{array}\right. } \end{aligned}$$

(6.1)

We look for a candidate solution to HJB in the form

$$\begin{aligned} {\tilde{w}}(t,x,y,d)= & {} A(T-t) (d-x)^2 + B(T-t) y^2 + F(T-t) (d-x)y \nonumber \\&+ \, G(T-t) (d-x) + H(T-t) y + K(T-t), \end{aligned}$$

(6.2)

for some deterministic functions A, B, F, G, H and K. Plugging the candidate function $\tilde{w}$ into equation (6.1), we see that ${\tilde{w}}$ is solution to the HJB equation iff the following system of ordinary differential equations (ODEs) is satisfied by A, B, F, G, H and K:

$$\begin{aligned} \left\{ \begin{array}{r} \textstyle A'+\frac{1}{4\gamma }(-2A+\nu F)^2 = 0 \\ \textstyle B'+\frac{1}{4\gamma }(2\nu B-F+1)^2 = 0 \\ \textstyle F'+\frac{1}{2\gamma } (-2A+\nu F) (2\nu B-F+1) = 0 \\ \textstyle G'-2\mu A +\frac{1}{2\gamma }(-2A+\nu F)(-G+\nu H) = 0 \\ \textstyle H'-\mu F +\frac{1}{2\gamma }(2\nu B-F+1)(-G+\nu H) = 0 \\ \textstyle K'-\mu G -(\sigma _0^2B+\sigma _d^2 A + \rho \sigma _0\sigma _d F)+\frac{1}{4\gamma }(-G+\nu H)^2 = 0 \end{array} \right. \end{aligned}$$

with the initial conditions $A(0)=\frac{1}{2}r(\eta , \beta )$, $B(0)=0$, $F(0)=0$, $G(0)=0$, $H(0)=0$, $K(0)=0$. We first solve the Riccati system relative to the triple (A, B, F), and obtain:

$$\begin{aligned} {\left\{ \begin{array}{ll} \textstyle A(t) = \frac{r(\eta ,\beta )(\frac{\nu }{2}t+\gamma )}{(r(\eta ,\beta )+\nu )t+2\gamma } , \\ \textstyle B(t) = -\frac{1}{2}\frac{t}{(r(\eta ,\beta )+\nu )t+2\gamma }, \quad F(t) = \frac{r(\eta ,\beta )t}{(r(\eta ,\beta )+\nu )t+2\gamma }. \end{array}\right. } \end{aligned}$$

(6.3)

Then we solve the first-order linear system of ODE relative to the pair (G, H), which leads to the explicit solution:

$$\begin{aligned} G(t) = 2\mu t A(t), \;\;&\hbox { and }&\; H(t) = -2 r(\eta ,\beta )\mu t B(t). \end{aligned}$$

(6.4)

Finally, we explicitly obtain K from the last equation:

$$\begin{aligned} K(t)= & {} \gamma \frac{\sigma _0^2 + \sigma _d^2r^2(\eta ,\beta ) - 2\rho \sigma _0\sigma _dr(\eta ,\beta )}{\big (r(\eta ,\beta ) + \nu \big )^2} \ln \Big ( 1 + \frac{(r(\eta ,\beta ) + \nu )t}{2\gamma }\Big ) \nonumber \\&+ \, \frac{\sigma _d^2r(\eta ,\beta )\nu +2\rho \sigma _0\sigma _dr(\eta ,\beta ) - \sigma _0^2}{2\big (r(\eta ,\beta ) + \nu \big )}t + \frac{r(\eta ,\beta )\mu ^2 t^2(\frac{\nu }{2}t+\gamma )}{ (r(\eta ,\beta ) + \nu )t+2\gamma }. \end{aligned}$$

(6.5)

By construction, ${\tilde{w}}$ in (6.2) with A, B, F, G, H and K explicitly given by (6.3)–(6.4)–(6.5), is a smooth solution with quadratic growth condition to the HJB equation (6.1). Moreover, the argmin in HJB equation for ${\tilde{w}}$ is attained for

$$\begin{aligned} {\tilde{q}}(t,x,y,d)= & {} - \frac{1}{2\gamma }\left[ \displaystyle \frac{\partial {\tilde{w}}}{\partial x} + \nu \displaystyle \frac{\partial {\tilde{w}}}{\partial y} + y \right] \\= & {} \frac{r(\eta ,\beta )(\mu (T-t)+d-x)-y}{(r(\eta ,\beta )+\nu )(T-t)+2\gamma } =: {\hat{q}}(T-t,d-x,y). \end{aligned}$$

Notice that ${\hat{q}}$ is linear, and Lipschitz in x, y, d, uniformly in time t, and so given an initial state (x, y, d) at time t, there exists a unique solution $({\hat{X}}^{t,x,y,d},{\hat{Y}}^{t,x,y,d},D^{t,d})_{t\le s\le T}$ to (2.1)–(2.7)–(2.8) with the feedback control ${\hat{q}}_s = {\hat{q}}(T-s,D_s^{t,d}-{\hat{X}}_s^{t,x,y,d},{\hat{Y}}_s^{t,x,y,d})$, which satisfies: $\mathbb {E}[\sup _{t\le s\le T} |{\hat{X}}_s^{t,x,y,d}|^2 + |{\hat{Y}}_s^{t,x,y,d}|^2 + |D_s^{t,d}|^2] < \infty$. This implies in particular that $\mathbb {E}[\int _t^T |{\hat{q}}_s|^2 ds] < \infty$, hence ${\hat{q}} \in \mathcal {A}_t$. We now call on a classical verification theorem (see e.g. Theorem 3.5.2 in [9]), which shows that ${\tilde{w}}$ is indeed equal to the value function ${\tilde{v}}$, and ${\hat{q}}$ is an optimal control. Finally, once the optimal trading rate ${\hat{q}}$ is determined, the optimal production is obtained from the optimization over $\xi \in \mathbb {R}$ of the terminal cost $C(D_T^{t,d}-{\hat{X}}_T^{t,x,y,d},\xi )$, hence given by: ${\hat{\xi }}_T = \frac{\eta }{\eta +\beta }(D_T^{t,d}-{\hat{X}}_T^{t,x,y,d})$. $\square$

1.2 Proof of Theorem 4.1

The Hamilton-Jacobi-Bellman (HJB) integro-differential equation arising from the dynamic programming associated to the stochastic control problem ${\tilde{v}} = {\tilde{v}}^{(\lambda )}$ with jumps in the dynamics of Y and D is:

$$\begin{aligned} {\left\{ \begin{array}{ll} \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial t} + \inf _{q \in \mathbb {R}} \Bigg [ q \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial x} + \nu q \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial y} + \mu \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial d} + \frac{1}{2} \sigma _0^2 \displaystyle \frac{\partial ^2 {\tilde{v}}^{(\lambda )}}{\partial y^2}\\ + \frac{1}{2} \sigma _d^2 \displaystyle \frac{\partial ^2 {\tilde{v}}^{(\lambda )}}{\partial d^2} + \rho \sigma _0\sigma _d \displaystyle \frac{\partial ^2 {\tilde{v}}^{(\lambda )}}{\partial y \partial d} + q(y+\gamma q) \Bigg ]\nonumber \\ + \lambda \big [ p^+{\tilde{v}}^{(\lambda )} (t,x,y+\pi ^+,d+\delta ^+) + p^- {\tilde{v}}^{(\lambda )}(t,x,y+\pi ^-,d +\delta ^-) - {\tilde{v}}^{(\lambda )}(t,x,y,d) \big ] = 0 \nonumber \\ {\tilde{v}}^{(\lambda )}(T,x,y,d) = {\tilde{C}}(d-x) = \frac{1}{2} r(\eta ,\beta ) (d-x)^2. \nonumber \end{array}\right. } \end{aligned}$$

Notice that with respect to the no jump case, there is in addition a linear integro-differential term in the HJB equation (which does not depend on the control), and the argmin is attained as in the no jump case for

$$\begin{aligned} {\tilde{q}}^{(\lambda )}(t,x,y,d)= & {} - \frac{1}{2\gamma }\left[ \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial x} + \nu \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial y} + y \right] . \end{aligned}$$

The HJB equation is then rewritten as

$$\begin{aligned} \!{\left\{ \begin{array}{ll}\! \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial t} + \mu \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial d} + \frac{1}{2} \sigma _0^2 \displaystyle \frac{\partial ^2 {\tilde{v}}^{(\lambda )}}{\partial y^2} + \frac{1}{2} \sigma _d^2 \displaystyle \frac{\partial ^2 {\tilde{v}}^{(\lambda )}}{\partial d^2} + \rho \sigma _0\sigma _d \displaystyle \frac{\partial ^2 {\tilde{v}}^{(\lambda )}}{\partial y \partial d} - \frac{1}{4\gamma }\left[ \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial x} + \nu \displaystyle \frac{\partial {\tilde{v}}^{(\lambda )}}{\partial y} + y \right] ^2 \\ + \lambda \big [ p^+{\tilde{v}}^{(\lambda )}(t,x,y+\pi ^+,d +\delta ^+) + p^- {\tilde{v}}^{(\lambda )}(t,x,y +\pi ^-,d+\delta ^-) - {\tilde{v}}^{(\lambda )}(t,x,y,d) \big ] = 0 \\ {\tilde{v}}^{(\lambda )}(T,x,y,d) = \frac{1}{2} r(\eta ,\beta ) (d-x)^2. \end{array}\right. } \end{aligned}$$

(6.6)

We look again for a candidate solution to (6.6) in the form

$$\begin{aligned} {\tilde{w}}^{(\lambda )}(t,x,y,d)= & {} A_{\lambda }(T-t) (d-x)^2 + B_{\lambda }(T-t) y^2 + F_{\lambda }(T-t) (d-x)y \nonumber \\&+ \, G_{\lambda }(T-t) (d-x) + H_{\lambda }(T-t) y + K_{\lambda }(T-t), \end{aligned}$$

(6.7)

for some deterministic functions $A_{\lambda }$, $B_{\lambda }$, $F_{\lambda }$, $G_{\lambda }$, $H_{\lambda }$ and $K_{\lambda }$. Plugging the candidate function ${\tilde{w}}^{(\lambda )}$ into equation (6.6), we see that ${\tilde{w}}^{(\lambda )}$ is solution to the HJB equation iff the following system of ordinary differential equations (ODEs) is satisfied by $A_{\lambda }$, $B_{\lambda }$, $F_{\lambda }$, $G_{\lambda }$, $H_{\lambda }$ and $K_{\lambda }$:

$$\begin{aligned} \left\{ \begin{array}{r} A_{\lambda }'+\frac{1}{4\gamma }(-2A_{\lambda }+\nu F_{\lambda })^2 = 0 \\ B_{\lambda }'+\frac{1}{4\gamma }(2\nu B_{\lambda }-F_{\lambda }+1)^2 = 0 \\ F_{\lambda }'+\frac{1}{2\gamma } (-2A_{\lambda }+\nu F_{\lambda }) (2\nu B_{\lambda }-F_{\lambda }+1) = 0 \\ G_{\lambda }'-2\mu A_{\lambda } +\frac{1}{2\gamma }(-2A_{\lambda }+\nu F_{\lambda })(-G_{\lambda }+\nu H_{\lambda })-\lambda (2\delta A_{\lambda }+\pi F_{\lambda }) = 0 \\ H_{\lambda }'-\mu F_{\lambda } +\frac{1}{2\gamma }(2\nu B_{\lambda }-F_{\lambda }+1)(-G_{\lambda }+\nu H_{\lambda })-\lambda (2\pi B_{\lambda } +\delta F_{\lambda }) = 0 \\ K_{\lambda }'-\mu G_{\lambda } -(\sigma _0^2 B_{\lambda } + \sigma _d^2 A_{\lambda } + \rho \sigma _0\sigma _d F_{\lambda })+\frac{1}{4\gamma }(-G_{\lambda } + \nu H_{\lambda })^2 \qquad \quad \\ \quad -\lambda [(p^+(\delta ^+)^2+p^-(\delta ^-)^2) A_{\lambda } + (p^+(\pi ^+)^2+p^-(\pi ^-)^2) B_{\lambda }\qquad \quad \\ +(p^+\delta ^+\pi ^++p^-\delta ^-\pi ^-) F_{\lambda } + \delta G_{\lambda } + \pi H_{\lambda }] = 0 \end{array} \right. \end{aligned}$$

with the initial conditions $A_{\lambda }(0)=\frac{1}{2}r(\eta , \beta )$, $B_{\lambda }(0)=0$, $F_{\lambda }(0)=0$, $G_{\lambda }(0)=0$, $H_{\lambda }(0)=0$, $K_{\lambda }(0)=0$. We first solve the Riccati system relative to the triple $(A_{\lambda },B_{\lambda },F_{\lambda })$, which is the same as in the no jump case, and therefore obtain: $A_{\lambda } = A$, $B_{\lambda } = B$, $F_{\lambda } = F$ as in (6.3). Then we solve the first-order linear system of ODE relative to the pair $(G_{\lambda },H_{\lambda })$, which involves the jump parameters $\lambda$, $\pi$ and $\delta$, and get:

$$\begin{aligned} G_{\lambda }(t)= & {} G(t) +\frac{\lambda }{2}\frac{r(\eta ,\beta )t(\pi t +2\delta (\nu t+2\gamma ))}{(r(\eta ,\beta )+\nu )t+2\gamma }, \\ H_{\lambda } (t)= & {} H(t) - \frac{\lambda }{2}\frac{(\pi - 2r(\eta ,\beta )\delta )t^2}{(r(\eta ,\beta )+\nu )t+2\gamma }, \end{aligned}$$

where G and H are given from the no jump case (6.4). Finally, after some tedious but straightforward calculations, we explicitly obtain $K_{\lambda }$ from the last equation:

$$\begin{aligned} K_{\lambda }(t)= & {} K(t)+ \lambda \gamma \frac{p^+(\pi ^+ - r(\eta ,\beta )\delta ^+)^2 +p^-(\pi ^- - r(\eta ,\beta )\delta ^-)^2}{\big (r(\eta ,\beta ) + \nu \big )^2} \ln \Big ( 1 + \frac{(r(\eta ,\beta ) + \nu )t}{2\gamma }\Big ) \\&- \, \frac{\lambda }{2} \frac{p^+((\pi ^+)^2 - r(\eta ,\beta ) \delta ^+(2\pi ^++\nu \delta ^+)) +p^-((\pi ^-)^2- r(\eta ,\beta ) \delta ^-(2\pi ^- +\nu \delta ^-))}{r(\eta ,\beta ) + \nu } t \\&+ \, \frac{\lambda r(\eta ,\beta )}{2}\frac{2\nu \mu \delta +\lambda ((p^+)^2\delta ^+(\pi ^++\nu \delta ^+)+(p^-)^2\delta ^-(\pi ^-+\nu \delta ^-))}{ r(\eta ,\beta ) + \nu }t^2 \\&+ \, \lambda ^2 \gamma r(\eta ,\beta ) \frac{r(\eta ,\beta )\delta ^2 +2\nu p^+p^-\delta ^+\delta ^- -((p^+)^2\delta ^+\pi ^++(p^-)^2\delta ^-\pi ^-)}{(r(\eta ,\beta )+ \nu )\big ((r(\eta ,\beta )+\nu )t + 2 \gamma \big )}t^2 \\&+ \, \frac{2\lambda \gamma r(\eta ,\beta )^2 \mu \delta }{(r(\eta ,\beta )+ \nu )\big ((r(\eta ,\beta )+\nu )t + 2 \gamma \big )}t^2 - \frac{\lambda ^2\pi ^2}{48\gamma }t^3 \\&+ \, \frac{\lambda ^2p^+p^-r(\eta ,\beta )}{2}\frac{2\nu \delta ^+\delta ^-+\delta ^-\pi ^++\delta ^+\pi ^-}{(r(\eta ,\beta )+\nu )t + 2\gamma }t^3\\&+ \, \frac{1}{8}\frac{4r(\eta ,\beta ) \mu \lambda \pi - \lambda ^2\pi ^2}{(r(\eta ,\beta )+\nu )t + 2\gamma }t^3, \end{aligned}$$

with K in (6.5). The function ${\tilde{w}}^{(\lambda )}$ in (6.7) may thus be rewritten as the sum of ${\tilde{w}}$ in (6.2) and another function of t, $d-x$ and y, and is by construction a smooth solution with quadratic growth condition to the HJB equation (6.6). Moreover, the argmin in HJB equation for ${\tilde{w}}^{(\lambda )}$ is attained for

$$\begin{aligned} {\tilde{q}}^{(\lambda )}(t,x,y,d)= & {} - \frac{1}{2\gamma }\left[ \displaystyle \frac{\partial {\tilde{w}}^{(\lambda )}}{\partial x} + \nu \displaystyle \frac{\partial {\tilde{w}}^{(\lambda )}}{\partial y} + y \right] \\= & {} \frac{r(\eta ,\beta )(\mu (T-t)+d-x)-y}{(r(\eta ,\beta )+\nu )(T-t)+2\gamma } \\&+ \, \lambda \frac{r(\eta ,\beta )\delta (T-t) + \frac{\pi }{4\gamma }(r(\eta ,\beta )+ \nu )(T-t)^2}{(r(\eta ,\beta )+\nu )(T-t) + 2\gamma } \\=: & {} {\hat{q}}^{(\lambda )}(T-t,d-x,y). \end{aligned}$$

Again, notice that ${\hat{q}}^{(\lambda )}$ is linear, and Lipschitz in x, y, d, uniformly in time t, and so given an initial state (x, y, d) at time t, there exists a unique solution $({\hat{X}}^{t,x,y,d},{\hat{Y}}^{t,x,y,d},D^{t,d})_{t\le s\le T}$ to (2.1)–(4.3)–(4.1) with the feedback control ${\hat{q}}_s^{(\lambda )} = {\hat{q}}^{(\lambda )}(T-s,D_s^{t,d}-{\hat{X}}_s^{t,x,y,d},{\hat{Y}}_s^{t,x,y,d})$, which satisfies: $\mathbb {E}[\sup _{t\le s\le T} |{\hat{X}}_s^{t,x,y,d}|^2 + |{\hat{Y}}_s^{t,x,y,d}|^2 + |D_s^{t,d}|^2] < \infty$, see e.g. Theorem 1.19 in [8]. This implies that $\mathbb {E}[\int _t^T |{\hat{q}}_s^{(\lambda )}|^2 ds] < \infty$, hence ${\hat{q}}^{(\lambda )} \in \mathcal {A}_t$. We now call on a classical verification theorem for stochastic control of jump-diffusion processes (see e.g. Theorem 3.1 in [8]), which shows that ${\tilde{w}}^{(\lambda )}$ is indeed equal to the value function ${\tilde{v}}^{(\lambda )}$, and ${\hat{q}}^{(\lambda )}$ is an optimal control. Finally, once the optimal trading rate ${\hat{q}}^{(\lambda )}$ is determined, the optimal production is obtained from the optimization over $\xi \in \mathbb {R}$ of the terminal cost $C(D_T^{t,d}-{\hat{X}}_T^{t,x,y,d},\xi )$, hence given by: ${\hat{\xi }}_T^{(\lambda )} = \frac{\eta }{\eta +\beta }(D_T^{t,d}-{\hat{X}}_T^{t,x,y,d})$. $\square$