Month: June 2011

VPP Pricing II: Mixed Integer Linear Programming

hpcquantlibLeave a comment

The next two steps are defining a simple VPP contract (or a simplified gas-run power plant) and setting up a mixed integer linear programming optimization (MIP) to calculated the intrinsic value and an upper bound for the extrinsic value based on a Monte-Carlo simulation and assuming perfect foresight. The third step outlined in the next part will then be the “exact” pricing of the extrinsic value using dynamic programming and finite difference methods.

The set-up of the simplified gas-run power plant is similar to the one explained in chapter 4.2.3 of the text-book [1]. In general the power plant has three power output level:

  • Plant is off, P=0
  • Generation at minimum load P=P_{min}
  • Generation at maximum load P=P_{max}

The power plant has a fixed efficiency rate

\zeta=\frac{MWh Power Output}{MWh Heat}.

Ramp rates will be neglected, but the power plant has a minimum uptime t_{up} and a minimum downtime t_{down}. The start-up costs are given by a fixed start-up cost \eta (in €) and the price of the gas needed to produce the start-up heat \theta (in MWh).

The mixed integer linear optimization is running in one hour blocks and is using three decision variables per hour i. The binary decision variable \beta_i is true if the power plant is running at minimum load P_{min} or at maximum load P_{max} and \beta_i is false if the plant is off. The real decision variable 0\le s_i \le 1 is equal to one if the plant is started in hour i, which is implied by the following constraint

\beta_i - \beta_{i-1} \le s_i.

The minimum up-time t_{up} and the minimum down-time t_{down} is a consequence of the constraints

\begin{array}{rcl} \beta_i &\ge& \sum_{t=i-t_{up}+1}^{t=i} s_t \\[7pt] \beta_i &\le& 1-\sum_{t=i+1}^{t=i+t_{down}} s_t \end{array}

The real decision variable 0\le \gamma \le 1 is equal to one if the power plant is running at maximum load P_{max} and zero if the power plant is either running at minimum load P_{min} or if the plant is off, that means

\beta_i \ge \gamma_i

Let P_i be the power price, G_i be the gas price and CO_2^i be the carbon dioxide price at hour i. The objective function is then given by

P\& L = \sum_{t=1}^N\left[\left(\gamma_iP_{max} + P_{min}(\beta_i-\gamma_i)\right) \left(P_i - \frac{G_i+CO_2^i}{\zeta}\right) - s_i\left(\eta + \theta (G_i+CO_2^i)\right)\right]

For a one year span the problem consists of 3\cdot 365 \cdot 24 = 26280 decision variables \{\beta_i, \gamma_i, s_i\} and 4 \cdot 365\cdot 24 = 35040 constraints. This comparable small problem can be solved using e.g. the Gnu Linear Programming Kit (GLPK). For an overview on open source linear/mixed integer programming solver see [2].

The model parameters and the example forward curves are outlined in the previous entry VPP Pricing I. The diagram below shows the intrinsic value and the upper bound for the total value (intrinsic plus extrinsic value) based on Monte-Carlo, perfect foresight and MIP for different  power plant efficiencies \zeta. The parameters of the VPP contract are given by

t_{up}=t_{down}=2h, P_{min}=8MW, P_{max}=40MW, \eta=300 EUR, \theta=20MWh,

the (fixed) carbon dioxide price is 3.0€ per MWh heat.

The source code is available here. It depends on GLPK and the latest QuantLib version from the SVN trunk or the next QuantLib 1.2 release.

It is now quite easy to add and price time-integral constraints, e.g. the following constraint restricts the number of starts within a year to be less than or equal to a given number

\sum_{t=1}^N s_i \le \#Starts.

The following diagram shows the results for \#Starts \le 25 and a minimum load P_{min}=25MW.

The source code is available here. It depends on QuantLib 1.1 and if  you want to generate the plot directly from the C++ program you’ll also need R, RCPP and RInside.

[1] M. Burger, B. Graeber, G. Schindlmayr, Managing Energy Risk, ISDN 978-0-470-ß2962-6

[2] S. R. Thorncraft, Evaluation of Open-Source LP Optimization Codes in Solving Electricity Spot Market Optimization Problems.

Advertisement

VPP Pricing I: Stochastic Processes & Partial Integro Differential Equation

hpcquantlib5 Comments

The basis for the virtual power plant (VPP) pricing algorithm is a joint stochastic process for the power price and the gas price. Instead of defining a  stochastic differential equation for the spark spread directly the spark spread is then given by the difference between the power price and the gas price times the VPP heating rate.

The Kluge model will be used to define the power price process [1], an exponential Ornstein-Uhlenbeck will describe the gas price [2]

\begin{array}{rcl} P_t &=& \exp(p_t + X_t + Y_t) \\ dX_t &=& -\alpha X_tdt + \sigma_x dW_t^x \\ dY_t &=& -\beta Y_{t-}dt+J_tdN_t \\ \omega(J) &=& \eta e^{-\eta J} \\ G_t &=& \exp(g_t + U_t) \\ dU_t &=& -\kappa U_tdt+\sigma_udW_t^u \\ \rho &=& \mathrm{corr} (dW_t^x, dW_t^u) \end{array}

where N_t is a Poisson process with jump intensity \lambda. To match the power forward curve F_0^t the seasonal function p_t is given by [1]:

p_t = \ln F_0^t -X_0 e^{-\alpha t}-Y_o e^{-\beta t} -\frac{\sigma_x^2}{4\alpha}\left(1-e^{-2\alpha t} \right ) - \frac{\lambda}{\beta}\ln\left( \frac{\eta-e^{-\beta t}}{\eta-1}\right), \eta\ge 1

To be consistent with the gas forward curve H_0^t the seasonal function g_t is defined by

g_t = \ln H_0^t -U_0 e^{-\kappa t}-\frac{\sigma_u^2}{4\kappa}\left(1-e^{-2\kappa t} \right )

The Feynman-Kac theorem can be applied to derive the corresponding three-dimensional  partial integro differential equation

\begin{array}{rcl} rV&=&\frac{\partial V}{\partial t}+\frac{\sigma_x^2}{2}\frac{\partial^2 V}{\partial x^2}-\alpha x\frac{\partial V}{\partial x}-\beta y \frac{\partial V}{\partial y} \\[6pt]&+&\frac{\sigma_u^2}{2}\frac{\partial^2 V}{\partial u^2}- \kappa u\frac{\partial V}{\partial u} +\rho\sigma_x\sigma_u\frac{\partial^2 V}{\partial x\partial u}\\[6pt] &+&\lambda\int_\mathbb{R}\left(V(x,y+z,u,t)-V(x,y,u,t) \right )\omega(z)dz \\ \end{array}

In general at least one further dimension is needed to keep track of the state of the virtual power plant. Therefore solving this model using finite difference methods will lead to a four-dimensional PIDE problem.

The following diagram shows an one year example path for the power and the gas price based on the freely available Kyos example forward curves. The sample parameterization is affected by [1] for the power process and [2] for the gas process.

\beta=200, \eta=2.5, \lambda=4, \alpha=7, \sigma_x=1.4, \kappa=4.45, \sigma_u=\sqrt{1.3}, \rho=70\%

The code is available here. It depends on the latest QuantLib version from the SVN trunk. If you want to generate the plot directly out of the C++ program you also need R, RCPP and RInside.

[1] T. Kluge, Pricing Swing Options and other Electricity Derivatives

[2] G. Fusai, A. Roncoroni, Implementing Models in Quantitative Finance: Models and Cases, Chapter 19, ISBN: 978-3-540-22348-1

Swing Option: Linear vs. Dynamic Programming II

hpcquantlib4 Comments

In order to get an improved impression on the differences between Monte-Carlo/linear programming and  finite difference methods/dynamic programming  a more realistic swing option with hourly payoff profile is evaluated based on a German hourly forward curve. The forward curve is taken from the Kyos example download page. The parameterization of the Kluge model is outlined in [1].

With maturity of twelve weeks the example swing option provides 12*7*24=2016 exercise opportunities. The number of overall exercised swing opportunities is constraint by

100 \le \sum_{i=1}^{2016} \beta_i \le 500.

The size of two dimensions of the finite difference method (dynamic programming) is therefore already been given. 2016 steps are needed in time direction and 500 steps are needed in the “consumed exercises” direction. Together with the two other directions – one for the power price and one dimension for the jump process – and without further simplidsfications this forms a pretty large finite difference problem. The Monte-Carlo based linear programming approach reduces the computational burden but will lead to an upper bound of the swing option price. The diagram below shows the corresponding results.



The code is available here. It depends on the GNU Linear Programming Kit, the Boost Thread library for parallelization and at the time of writing on the latest QuantLib version from the SVN trunk. If you want to generate the plot directly out of the C++ program you also need R, RCPP and RInside. To utilize all CPU cores please use the -DNTHREADS=(number of CPU cores)  compiler switch. In addition to run the program you have to download the forward curve German power from the Kyos web page and convert it into a text file EEX_2010-2013.txt of the folliwing format

4-Oct-2010 36.73 32.09 27.32 23.22 25.71 35.60 47.32 …

5-Oct-2010 38.12 31.12 22.76 25.65 27.87 34.60 50.01 …

[1] T. Kluge, Pricing Swing Options and other Electricity Derivatives