Numerical Methods for Valuation and Optimal Operation of Natural Gas Storage

(1)

Numerical Methods for Valuation and Optimal Operation of Natural Gas Storage

Erik Magnus G. Følstad

Master of Science in Physics and Mathematics Supervisor: Espen Robstad Jakobsen, MATH

Department of Mathematical Sciences Submission date: March 2015

Norwegian University of Science and Technology

(2)

(3)

Abstract

The thesis describes different approaches for solving numerically a PDE model for the valuation and optimal operation of natural gas storage, characterized as a Hamilton Jacobi Bellman (HJB) equation. The HJB equation is derived by for- mulating the given natural gas storage problem as a stochastic control problem and then applying the dynamic programming principle. We present three separate numerical methods for solving the HJB equation, namely a standard upwind finite difference method, and two new methods characterized as: (i) a semi-Lagrangian time stepping method combined with a one dimensional finite element method, and (ii) a two dimensional finite element method combined with finite difference discretization in time. The upwind finite difference method is shown to be consistent, stable and monotone. These properties guarantee that the numerical solution converge to the viscosity solution of the HJB equation, [19]. Numerical results suggest that the two new methods converge to the same solution as the finite difference method for a given test case.

(4)

Sammendrag

Denne oppgaven beskriver ulike numeriske metoder for ˚a løse en PDE modell for verdisetting og optimal styring av naturgasslagring, gitt som en Hamilton Jacobi Bellman (HJB) likning. Likningen er utledet ved ˚a først formulere problemstillin- gen som et stokastisk kontrollproblem og deretter benytte dynamisk programmer- ingsprinsippet. Vi presenterer tre ulike metoder for ˚a løse denne likningen, nærmere bestemt en standard oppvind differanse metode, og to nye metoder beskrevet som:

(i) en semi-Lagrangian tids-diskretiseringsmetode kombinert med endelig element metode i en romlig retning, og (ii) en endelig elementmetode i to romlige retninger kombinert med endelig differanse i tid. Det blir vist at differansemetoden er konsis- tent, stabil og monoton, hvilket impliserer at den numeriske løsningen konvergerer til viskositetsløsningen av HJB likningen, [19]. Numeriske resultater tyder p˚a at de to nye metodene konvergerer mot samme løsning som differansemetoden for et gitt test-problem.

(5)

(6)

List of Figures

4.1 Node placement on a generic elementKi= [xi, xi+1]∈ Th, for linear, quadratic and cubic elements. Nodes displayed as black dots. . . 35 6.1 P1, P2 and P3 triangular elements with nodes (degrees of freedom)

displayed in black. . . 50 6.2 Mapping ˆF(ˆx,y) := (g(ˆˆ x,y), h(ˆˆ x,y))ˆ ^> from the reference element ˆK

(left) to a physical elementK (right). . . 54 7.1 Numerical solution of the value function and the optimal control

policy evaluated atτ = 1. The discretization parameters are ∆x= xmax/20,∆y=ymax/20 and ∆τ=T /200. . . 65 7.2 Numerical solution of the value function after 10 time steps, with

h= 0.04 and ∆τ = 0.0156. When we remove the edge stabilization, spurious oscillations are observed in the solution. The oscillations starts inx= 0, where there is least diffusion. Similar observation are made withP1elements. We found that the stabilization parameters γ= 0.02 andγ= 0.5 worked well forP1 andP2 elements respectively. 66

(9)

(10)

List of Tables

4.1 Gauss quadrature rules for the reference interval (−1,1). . . 37 6.1 Gauss quadrature rules for the reference triangle. . . 54 7.1 Verification of the convergence rate ˆβ = 2 for the p1-FE method.

withα= 0.7071. The stabilization parameter isγ= 0.5 . . . 62 7.2 Verification of the convergence rate ˆβ = 3 for the p2-FE method.

withα= 0.7937. The stabilization parameter isγ= 0.02 . . . 62 7.3 Verification of the convergence rateβ = 2 for thep1-SLFE method

withα= 0.7071. No edge stabilization is used (γ= 0). . . 63 7.4 Verification of the convergence rateβ = 3 for thep2-SLFE method

withα= 0.7937. No edge stabilization is used (γ= 0). . . 63 7.5 Verification of convergence to the viscosity solution for theP1-SLFE

method; ˜vk and ˆvk represent the solutions of theP1-SLFE-method and the finite difference method, respectively. No edge stabilization is used (γ= 0). . . 64 7.6 Verification of convergence to the viscosity solution for theP2-SLFE

method; ˜vk and ˆvk represent the solutions of theP2-SLFE-method and the finite difference method, respectively. No edge stabilization is used (γ= 0). . . 64 7.7 Verification of convergence to the viscosity solution for the P1-FE

method; ˜v_k and ˆv_k represent the solutions of theP1-FE-method and the finite difference method, respectively. The stabilization parameter isγ= 0.5. . . 64 7.8 Verification of convergence to the viscosity solution for the P2-FE

method; ˜vk and ˆvkrepresent the solutions of the FE-method and the finite difference method, respectively. The stabilization parameter isγ= 0.02. . . 64

(11)

(12)

Chapter 1

Introduction

The real option of storing natural gas in the presence of uncertain gas prices is important in planning-models for the production of natural gas and allows industries to exploit spot market variations and seasonal trends, [15]. As we will demonstrate, the real option of storing natural gas can be modeled as a stochastic control problem from which a Hamilton Jacobi Bellman (HJB) equation is derived. This thesis is concerned with the numerical solution of such models.

1.1 Previous Work

Thompson et al. [14] derive a natural gas storage model in the form of a HJB equation and propose a fully explicit finite difference scheme. Because of the hyperbolic nature of the HJB equation, the need for stabilization techniques to prevent numerical oscillations are addressed. Thompson et al. suggest a min-mod slope limiting method to cope with this problem.

Forsyth and Chen [19] present an implicit semi-Lagrange time stepping scheme combined with a finite difference method for the natural gas problem and show that the scheme converges to the viscosity solution of the HJB equation by proving that the scheme is consistent, stable and monotone. The semi-Lagrange time stepping method is widely used in the field of numerical weather predictions [13] and is known to be stabilizing for convection dominated problems.

1.2 The purpose of this thesis

• Provide modeling background and derive the HJB equation for the natural gas storage problem.

• Study numerical methods for solving the given HJB equation and propose new methods.

• Implement the methods in MATLAB and conduct numerical experiments.

(13)

1.3 Contribution

After providing a modeling background for the HJB equation, we study an upwind finite difference method which is shown to be consistent, stable and monotone provided a linear CFL-condition. Then two new methods for this problem are introduced:

(i) A semi-Lagrangian time stepping method in one spatial direction combined with a finite element method in the other spatial direction. The method is inspired from [19]. Another method that combines semi-Lagrangian time stepping with finite elements is given in [18], however, in [18] finite elements and semi-Lagrangian time stepping are applied in the same spatial direction.

(ii) A two dimensional finite element method combined with a finite difference discretization in time. To cope with the issue of spurious oscillations we have implemented the edge stabilization technique as given by Burman and Hansbo [4].

1.4 Outline

Chapter 2 Some basic ideas from the field of stochastic control theory are introduced. These ideas are applied in section 2.3 to derive a model for the valuation and optimal operation of a natural gas storage facility in the form of a Hamilton Jacobi Bellman equation.

Chapter 3 An upwind finite difference scheme is described and analyzed. The scheme is shown to be consistent, stable and monotone provided a linear CFL-condition.

Chapter 4 This chapter is to intended as a brief introduction to the finite element method which will be used in the chapter 5.

Chapter 5 A semi-Lagrange time stepping combined with a finite element method is presented. A fully implicit scheme is achieved via a linearization technique.

Chapter 6 We present a finite element method coupled with finite finite difference time stepping. Edge stabilization [4] is implemented to prevent numerical oscillations.

Chapter 7 Numerical experiments are conducted with a standard test case in the literature.

Chapter 8 We summarize the work and give some concluding remarks.

(14)

Chapter 2

Background

2.1 Stochastic control theory

Many real life phenomena, in which the random nature of certain quantities plays a significant role, can be successfully modeled as a stochastic differential equation (SDE). For instance, in finance, stochastic differential equations can be used to describe the evolution of stock prices or commodity prices [12].

In the field of Stochastic control theory, one studies how the state of a system, given as the solution of a system of SDE’s, is influenced by some controlled parameter. The objective is typically to maximize (or minimize) some functional that depends on the the state of the system by choosing the optimal value for the controlled parameter. As described in [9, p. 27], a stochastic control problem is characterized by the following features;

• State of the system: We consider a systemX= (X1, X2. . . , Xn) satisfying a system of SDE’s. The state of the system at any times≥0 is given byX(s).

• Control: The systemXis influenced by a controlled processα:s→A⊂R. The value ofαat timesmust be chosen by using only available information at times. In addition,αmust satisfy certain constraints to beadmissible.

• Performance criterion: The objective is to maximize (or minimize) some functional P(X, α) over all admissible controls. Hence, a stochastic control problem can be defined as finding an optimal controlα^∗, such thatP(X, α^∗) is maximized (or minimized).

Lett, T ∈R, such that 0≤t < T. Let f :Rⁿ×A→Rⁿ be a given function and letArepresent the set of admissible control functions

A={α: [t, T]→A|α(·) admissible}.

The dynamics of the state X of the system is dependent on the control process

(15)

α(·)∈ Aand is given by the SDE

dX(s) =f(X(s), α(s), s)ds+σ(X(s), s)dW(s), t < s < T (2.1)

X(t) =x, (2.2)

where dW(s) represents the standard increment of an n-dimensional Brownian motion and σis assumed to be a known deterministic function. The point (x, t) is referred to as the initial state of the system.

A generic performance functional can be stated as P[α(·)]x,t= E

"

Z T t

r(X(s), α(s), s)ds

X(t) =x

#

, ∀α∈ A. (2.3)

where ris some given function that typically represents cash flow or the instanta- neous rate of change in some other desirable quantity that accumulates over the time horizon [0, T]. The requirements for a functionα:T →Rto be admissible is problem dependent, e.g.Acan depend on the initial state of the system.

In this work, we will consider the following stochastic optimal control problem;

For all initial states (x, t)∈Ω×[0, T], findα^∗∈ A such that P[α^∗(·)]x,t= sup

α∈A

P[α(·)]x,t, (2.4)

with Ω ⊂ Rⁿ representing all possible values for x. The value function v : Ω× [0, T]→Ris as

v(x, t) = sup

α∈A

P[α^∗(·)]_x,t, ∀(x, t)∈Ω×[0, T]. (2.5) Under certain assumptions, one can show viadynamic programming that the value function satisfies the Hamilton Jacobi Bellman equation, which will be introduced in the next section.

2.2 Dynamic Programming Principle and the Hamilton Jacobi Bellman Equation

Bellman’s principle of optimality, also referred to as the dynamic programming principle, see [2, p. 83], is stated as

Dynamic Programming Principle. An optimal policy has the property that whatever the initial state and initial decision are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision.

With respect to the optimal control problem (2.4) defined in the previous section, the dynamic programming principle can be formulated mathematically as

v(x, t) = sup

α∈A

E

"

Z t+δt t

(r(X(s), α(s), s) +v(Xα(t+δt), t+δt))ds

#

, (2.6)

(16)

where X_α(t+δt) represents the state of the system X at time t+δt and the subscript α indicates that we have applied the control α(·) in the time interval [t, t+δt].

It can be shown that the value functionvsatisfies the Hamilton Jacobi Bellman equation, stated as

−∂v

∂t(x, t)−sup

a∈A

[L^av(x, t) +r(x, a, t)] = 0, (2.7) whereL^a is an operator including partial derivatives up to second order depending on the control parametera∈A. Equation (2.7) is derived from problem (2.4) using the dynamic programming principle (2.6) and assuming that the value function is sufficiently smooth [9]. We shall heuristically derive an instance of the Hamilton Jacobi Equation in section 2.3 relating to the optimal control of a natural gas storage facility.

2.3 Derivation of the Natural Gas Storage Model

We derive the Hamilton Bellman Jacobi equation for the valuation and optimal operation of a gas storage facility, following Thompson et al. [14]. Similar models can also be derived for the valuation and optimal operation of a hydro electrical power plant [8], [5].

For the purposes of this model we think of a natural gas storage facility as being essentially a storage unit from which gas can be injected and withdrawn at any time. In order to maximize profits the operators of the facility wish to, roughly speaking, sell gas when the price is high and buy gas when the price is low. We introduce the following variables

• X - the price per unit of gas,

• Y - the amount of working gas in storage.

It is assumed that the facility can only be operated within a finite time horizon.

We lett denote the present time and we letT denote the end of the time horizon.

For any time s ∈[t, T], the price of unit gas and the amount of gas in storage is expressed as X(s) andY(s), respectively.

The policy chosen by the operators of the facility is represented by a control function α: [t, T] →R. That is, for each times ∈[t, T] the value of the control functionα(s) represents the volume rate of gas being sold. We use the convention that ifα >0, then gas is withdrawn from the facility and sold. On the other hand ifα <0, then gas is being bought and injected into storage. For a control function α(·) to beadmissible we require that

amin(Y(s))≤α(s)≤amax(Y(s)), ∀s∈[t, T],

whereaminandamaxare known functions ofY representing the maximal injection rate and the maximal withdrawal rate respectively. In the following we let

A(y) ={a∈R:amin(y)≤a≤amax(y)}, ,

(17)

and we let A(y) represent the collection of all admissible control functions given the initial level of gasY(t) =y.

We assume that there is a leakage of gas represented by a deterministic function λ, perhaps dependent on the control policy α and the amount of gas in storage.

The flow of gas out from storage, represented by the functionf, is given by f(Y(s), α(s), s) =λ(Y(s), α(s), s) +α(s) (2.8) We are using the sign convention that iff >0 then gas is flowing out from storage and if f <0 then gas is flowing into storage. Consequently, the rate of change in the amount of gas in storage is given by ^dY_ds =−f. Lety represent the amount of working gas presently in storage at timet. ThenY satisfies the following ODE, (in infinitesimal form) :

dY(s) =−f(Y(s), α(s), s)ds t < s < T

Y(t) =y. (2.9)

The price per unit of gasX is assumed to be a stochastic process satisfying an SDE of the form

dX(s) =µ(X(s), s)ds+σ(X(s), s)dW(s), t < s < T

X(t) =x, (2.10)

where W represents aWiener process.

DEFINITION. A real valued stochastic processW(t)is called a Wiener process, or Brownian motion, if

1. W(0) = 0

2. each sample path is continuous

3. W(t)isN(0, t), (W is normally distributed with mean 0 and variancet) 4. W has independent increments

Therevenueof the gas storage facility is defined as the present value of the sum of allcash flow over the whole contract period. The cash flow, denotedr, is equal to the value of the gas currently being bought or sold minus the value of the gas currently being lost, that is

r(X(s), Y(s), α(s), s) = α(s)−λ(Y(s), α(s), s)

X(s), ∀s∈[t, T]. (2.11) Thepresent value of this cash flow is assumed to be given by

e^−ρ(s−t)(α(s)−λ(Y(s), α(s), s))X(s),

where ρ denotes therisk free interest rate [6, p. 115]. Consequently, the present value of the the sum of all cash flow over the whole contract period is equal to

Z T t

e^−ρ(s−t)r(X(s), Y(s), α(s), s)ds.

(18)

The performance functional is given as P[α(·)]x,y,t= E

"

Z T t

e^−ρ(s−t)r(X(s), Y(s), α(s), s)ds

X(t) =x, Y(t) =y

# ,

∀α∈ A(y), (2.12) which associates with each admissible control α ∈ A(y) an expected revenue P[α]x,y,t, given the initial stateX(t) =x, Y(t) =y.

Remark 2.1. Note that with respect to generic performance functional (2.3), the performance functional given by the previous equation has an extra discount factor e^−ρ(s−t).

We are interested in the following stochastic control problem for each initial state (x, y, t):

findα∈ A(y): P[α^∗(·)]x,y,t= sup

α∈A(y)

P[α(·)]x,y,t. (2.13) Recall that the value function is given as

v(x, y, t) = sup

α∈A(y)

P[α(·)]x,y,t, ∀(x, y, t)∈Ω×[0, T]. (2.14) We will now proceed to derive heuristically a partial differential equation for the value function v, following [9, p. 43], specifically we obtain the Hamilton Jacobi Bellman equation (2.7). We will be needing the following well known result from the theory of stochastic differential equations, as given in [6, p. 78].

ITO’S FORMULA. Let X= (X1, X2, . . . , Xn)represent an n-dimensional stochastic process such that

dXi=µi(X(s), s)ds+

n

X

j=1

Gij(X(s), s)dWj fori= 1, . . . , n.

with µ_i(X(s), s)∈L¹(0, T)andG_ij(X(s), s)∈L²(0, T).

Ifu:Rⁿ×[0, T]→R is continuous and the partial derivatives ^∂u_∂t,_∂x^∂u

i,_∂x^∂²^u

i∂xj, (i, j= 1, . . . , n) exist and are continuous, then

d(u(X₁, . . . , X_n, s)) = ∂u

∂tds+

n

X

i=1

∂u

∂x_idX_i+1 2

n

X

i,j=1

∂²u

∂x_i∂x_j

n

X

l=1

G_ilG_jlds.

Suppose that that the operators of the gas storage facility applies the constant control policy α≡a∈ A(y) within the time interval [t, t+δt]. According to the

(19)

dynamic programming principle (2.6), the value function satisfies the inequality v(x, y, t)≥E

"

Z t+δt t

e^−ρ(s−t)r(X, Y, a, s)ds+ e^−δtv(X(t+δt), Y(t+δt), t+δt) ds

# . (2.15) Provided that the value functionv is a sufficiently smooth function ofx, y andt, Ito’s formula implies that

dv(X(s), Y(s), s) = ∂v

∂t +µ∂v

∂x + 1 2σ²∂²v

∂x² −f ∂v

∂y

ds+σ∂v

∂xdW(s). (2.16) Remark 2.2. For notational convenience we have omitted to write out that the partial derivatives ^∂v_∂t,^∂v_∂x,^∂_∂x²^v2, ^∂v_∂y appearing on the right hand side of the last equation are evaluated at the point (X(s), Y(s), s), the functions µ,σ are evaluated at (X(s), s)andf is evaluated at(X(s), Y(s), a, s).

Inequality (2.15) can be rearranged as 0≥E

"

Z t+δt t

e^−ρ(s−t)r ds+ e^−ρδtδv−(1−e^−ρδt)v

#

. (2.17)

with

δv:=v(X(t+δt), Y(t+δt), t+δt)−v(x, y, t).

From equation (2.16) we obtain δv=

Z t+δt t

dv= Z t+δt

t

∂v

∂t +µ∂v

∂x +1 2σ²∂²v

∂x²−f∂v

∂y

ds+ Z t+δt

t

σ∂v

∂xdW(s).

By substituting the expression for δv given by the previous equation into (2.17) and using that

E

"

Z t+δt t

σ∂v

∂xdW(s)

#

= 0, we obtain

0≥E Z t+δt

t

e^−ρ(s−t)r ds + e^−ρδt

Z t+δt t

∂v

∂t +µ∂v

∂x+1 2σ²∂²v

∂x² −f∂v

∂y

ds−(1−e^−ρδt)v

. Divide the previous inequality byδtand take the limit asδtgoes to zero. Assuming that lim

δt→0and E[·] are interchangeable, we get 0≥E

δt→0lim 1 δt

Z t+δt t

e^−ρ(s−t)r ds + e^−ρδt

Z t+δt t

∂v

∂t +µ∂v

∂x +1 2σ²∂²v

∂x² −f∂v

∂y

ds−(1−e^−ρδt)v

.

(20)

Via the mean-value theorem, the previous equation implies that 0≥

∂v

∂t +1 2σ²∂²v

∂x² +µ∂v

∂x +ρv

(x, y, t) +r(x, y, a, t)−f(y, a, s)∂v

∂y(x, y, t).

(2.18) To obtain the last inequality we have also used the identity

lim

δt→0

1

δt(1−e^−ρδt) =ρ.

Since (2.18) is true for any admissible a∈A(y), we have 0≥

∂v

∂t +1 2σ²∂²v

∂x² +µ∂v

∂x −ρv

(x, y, t) + sup

a∈A(y)

r−f∂v

∂y

(x, y, a, t). (2.19) On the other hand, if the operators of the gas storage facility applies the optimal policy α^∗ within the interval [t, t+δt], it can be shown trough similar arguments that

0 = ∂v

∂t +1 2σ²∂²v

∂x² +µ∂v

∂x −ρv

(x, y, t) +r(x, y, a^∗, t)−f(y, a^∗, t)∂v

∂y(x, y, t), (2.20) with

a^∗:=α^∗(t).

Inequality (2.19) and equation (2.20) together imply that vshould satisfy 0 =

∂v

∂t +1 2σ²∂²v

∂x² +µ∂v

∂x−ρv

(x, y, t) + sup

a∈A(y)

r(x, y, a, t)−f(y, a, t)∂v

∂y(x, y, t)

.

The previous equation is formulated backwards in time, via the change of variable τ =T −t, we get

∂v

∂τ −1 2σ²∂²v

∂x² −µ∂v

∂x+ρv

(x, y, τ)

= sup

a∈A(y)

r(x, y, a, τ)−f(y, a, τ)∂v

∂y(x, y, t)

. (2.21)

2.4 Important assumptions about the NGS model

2.4.1 Boundary

In reality there is no upper bound on the price variablex. However, in our numerical solution of (2.21) we can only consider a finite price range, so let xmax represent

(21)

the maximal allowed price. The gas inventory variable y has an upper bound, represented by y_max, equal to the maximum storage capacity. In this work we consider solving equation (2.21) on the domain [0, x_max]×[0, y_max]×[0, T] and we refer to

Ω = [0, xmax]×[0, ymax],

as thespatial domain. For convenience we assume through out this work that the coefficients of (2.21) are scaled such that the spatial domain is equal to the unit square, that is

xmax=ymax= 1.

The gas price process is assumed to be mean reverting. Let µ0 > 0 and let

¯

x∈(0,1). A simple example of a mean reverting stochastic process is stated dX(s) =µ₀(X(s)−x)¯ ds+σ(X(s), s)dW(s),

which corresponds to setting

µ(X(s), s) =µ₀(X(s)−x),¯

in equation (2.10). The process is called mean reverting because the gas price is constantly drawn to the mean reversion level ¯x∈(0,1) which corresponds to the long term average market price. When the price process is mean reverting, we have the following conditions on the boundariesx= 0 andx= 1

µ|x=0≥0,

µ|x=1≤0. (2.22)

For the gas inventory process it is safe to assume that gas can only be pumped into storage if the current level of gas is less than the maximum capacity. On the other hand, gas can only be pumped out of the facility if there is any gas left in inventory. This reasoning leads to the conditions

f|y=0≤0

f|y=1≥0, (2.23)

where we recall thatf represents the flow of gas out from storage, (2.8). Negative values forf corresponds to gas being pumped in to storage and positive values for f corresponds to gas being released.

In [14] and [19] it is argued without any real justification that

∂²v

∂x² →0 as x→xmax. (2.24)

However, the boundary x = xmax := 1 is assumed to be “far away” from the realistic price range, so this condition might not have a big impact on the solution.

Condition (2.24) implies that the term σ²^∂_∂x²^v2 in equation (2.21) goes to zero as

(22)

x→1. As suggested by [5] this can be implemented by replacing σ with ˆσsuch that

ˆ

σ(x) =σ(x), ifx∈[0,1−],

ˆ

σ(x) = 0, ifx∈[1−,1],

where represents some small number. We will simply assume that σ → 0 as x→xmax. It is also assumed that σ→0 asx→0. Consequently, equation (2.21) is nearly hyperbolic close to the boundaries x= 0 and x= 1, (because the term σ²^∂_∂x²^v2 vanishes), and boundary conditions needs only be prescribed at the inflow part of the boundariesx= 0 andx= 1. However, condition (2.22) imply that the boundaries x= 0 and x= 1 are outflow boundaries. That is, the flow field given as

f = (µ,−f)^>,

always points inwards to the domain. Hence, no boundary conditions are needed for equation (2.21) at x ∈ {0,1}. Because the diffusive term in (2.21) has no component in the y-direction and the boundaries y = 0 and y = 1 are outflow boundaries by condition (2.23), no boundary conditions are needed at y∈ {0,1}.

For simplicity we will apply the initial condition v(x, y,0) = 0, as suggested in [14]. Some alternative initial conditions are discussed in [19].

2.4.2 Bang-bang controls

The model can be split into two equations as follows a^∗= arg sup

a∈A(y)

r(x, y, a^∗, τ)−f(y, a^∗, τ)∂v

∂y(x, y, t)

(2.25)

∂v

∂τ = 1

2σ²∂²v

∂x² +µ∂v

∂x−ρv

(x, y, τ)−f(y, a^∗, τ)∂v

∂y(x, y, t) +r(x, y, a^∗, τ) (2.26) In [19] and [14] the chosen model for the gas leakageλis of the form

λ(a) =

0 ifa≥0 k ifa <0 . With this choice ofλ, we get

r(x, a) = (a−λ(a))x, f(a) = (a+λ(a)),

(23)

so the feedback controla^∗can be found by solving two linear sub-problems, corre- sponding toa≥0 anda <0 respectively, given as

a^∗₁= arg sup

a∈A(y)

ax−a∂v

∂y(x, y, t)

, a≥0

a^∗₂= arg sup

a∈A(y)

(a−k)x−(a+k)∂v

∂y(x, y, t)

a <0

and then choose the best of these computed solutions:

a^∗= arg sup{(a^∗₁−λ(a^∗₁))x−(a^∗₁+λ(a^∗₁)),(a^∗₂−λ(a^∗₂))x−(a^∗₂+λ(a^∗₂))}

Since both sub-problems are linear inait is easy to verify that the optimal control only takes on values from the finite set{amin(y),0, amax(y)}. Hence, the optimiza- tion problem is solved by simply comparing the possible values of the expression r(x, a)−f(a)^∂v_∂y fora∈ {amin(y),0, amax(y)}. That is,

a^∗(x, y, τ) = arg sup

a∈{amin,0,a_max}

((a−λ(a))x−(a+λ(a)).

The controls are said to be of bang-bang type.

(24)

Chapter 3

A Monotone Finite Difference Scheme

We present a first order semi implicit finite difference scheme for solving equation (2.21). In the spatial discretization we employ the standard first order upwinding technique, as described in [10, p. 284]. The time discretization is implicit in the price direction (x-direction) and explicit in the inventory direction (y-direction), leading to a linear CFL-condition.

Even if the following method is only first order accurate it has other desirable properties. The scheme is l_∞-stable, consistent and monotone. As noted in [19], these properties ensures that the numerical scheme converges to theviscosity solu- tionof (2.21), which is the appropriate solution for stochastic control problems [9].

In addition, the finite difference method is generally easy to implement and vec- torization of the method in MATLAB is straight forward.

We refer to [16] for an upwind finite volume method for the Hamilton Jacobi Bellman equation.

3.1 Grid

Let (xi)^I_i=1, (yj)^J_j=1 and (τn)^N_n=0, be sets ofgrid points such that

0 =x1< x1<· · ·< xI−1< xI =xmax:= 1, (3.1) 0 =y1< y1<· · ·< yJ−1< yJ=ymax:= 1, (3.2) and

0 =τ0< τ1<· · ·< τN−1< τN =T, (3.3) Fro simplicity, we assume that the grid points are uniformly spaced, that is

∆x=xi+1−xi, i= 1, . . . , I,

(25)

∆y=y_j+1−y_j, j= 1, . . . , J, and

∆τ=τ_n+1−τ_n, n= 0, . . . , N.

for some ∆x >0, ∆y >0 and ∆τ >0.

3.2 The upwinding technique

In order to introduce the upwind method, consider the following linear convection diffusion equation

−σ∂²v

∂x² +b(x)∂v

∂x = 0. (3.4)

Let ˜v represent the numerical approximation ofv to be defined. According to the first order upwind method, the approximation ˜vxof the convective termvxis given by the following rule:

˜

v_x= ˜v(x)−˜v(x−∆x)

∆x ,ifb >0,

˜

vx= ˜v(x+ ∆x)−v(x)˜

∆x ,ifb <0.

(3.5)

Let

vi:= ˜v(xi), fori= 1, . . . , I, define the operators (·)⁺ and (·)⁻ such that

a⁺= max(a,0), a⁻= min(a,0), ∀a∈R, and let the difference operators D⁺_x and D⁻_x be given as

D⁺_x˜v(x) = 1

∆x(˜v(x+ ∆x)−v(x)),˜ D⁻_x˜v(x) = 1

∆x(˜v(x)−v(x˜ −∆x)).

By using the standard central difference approach to approximate the diffusive term in equation (3.4) and the upwind method, given by (3.5), to approximate the convective terms, we end up with the following numerical scheme:

−σD⁺_x(D⁻_xv_i) +b(x_i)⁺D⁻_xv_i+b(x_i)⁻D⁺_xv_i= 0, i= 1. . . , I−1.

Remark 3.1. We observe that

D⁺_x(D⁻_xvi) = vi+1−2vi+vi−1

∆x² .

(26)

3.3 Numerical scheme

In the previous section we introduced the upwind method for a one dimensional linear convection diffusion equation. We will now proceed to apply the method on equation (2.21), which we restate here for convenience:

∂v

∂τ −1 2σ²∂²v

∂x²−µ∂v

∂x +ρv

(x, y, τ) = sup

a∈A(y)

r(x, y, a, τ)−f(y, a, τ)∂v

∂y(x, y, t)

. (3.6) We assume, for simplicity, thatσandµonly depend onx. In what follows, we will use the notation

µ_i=µ(x_i), σ_i=σ(x_i), and

f_jⁿ⁺¹(a) =f(y_j, τ_n+1, a), r_ijⁿ⁺¹(a) =r(x_i, y_j, τ_n+1, a).

Let the difference operators D⁺_y and D⁻_y be defined such that

D⁺_yv(x, y, τ) =v(x, y+ ∆y, τ)−v(x, y, τ), D⁻_yv(x, y, τ) =v(x)−v(x, y−∆y, τ).

By following the same discretization procedure as in the previous section for each direction x,y andτ, we obtain from equation (3.6) the numerical scheme:

vⁿ⁺¹_i,j −v_i,jⁿ

∆τ −σ_iD⁺_x(D⁻_xvⁿ⁺¹_i,i )−µ⁺_iD⁺_xv_i,jⁿ⁺¹−µ⁻_i D⁻_xvⁿ⁺¹_i,j = max

a∈A(yj)

rⁿ_ij(a)−f_jⁿ(a)⁺D⁻v_i,jⁿ −f_jⁿ(a)⁻D⁺_yv_i,jⁿ , (3.7) fori= 0, . . . , I,j= 0, . . . , J andn= 0, . . . , N−1.

Remark 3.2. We have treated the convective term f^∂v_∂y and the term r(x, y, a, τ) in equation (3.6) explicitly, that is, in equation (3.7), these terms are evaluated at τ =τ_n instead of τ=τ_n+1.

Remark 3.3. As we will show in the next subsection, no boundary conditions are needed for (3.7), provided thatσ(0) =σ(1) = 0,µ(0)≥0,µ(1)≤0,f_|_y=0 ≤0and f_|_y=1≥0.

3.3.1 Boundary

At the boundaries x= 1 and x= 0, we apply the conditions

σ1= 0, σI = 0,

(27)

(recall section 2.4.1), and it is assumed that

µ₁≥0, µ_I ≤0,

so for i= 1 andi=I, respectively, the scheme given by (3.7) reads vⁿ⁺¹_1,j −vⁿ_1,j

∆τ −µ⁺₁D⁺_xv_1,jⁿ⁺¹= max

a∈A(yj) r_1,jⁿ (a)−f_jⁿ(a)⁺D⁻v_1,jⁿ −f_jⁿ(a)⁻D⁺_yvⁿ_1,j , (3.8) and

v_I,jⁿ⁺¹−v_I,jⁿ

∆τ −µ⁻₀D⁻_xv_I,jⁿ⁺¹= max

a∈A(yj)

r_I,jⁿ (a)−f_jⁿ(a)⁺D⁻v_I,jⁿ −f_jⁿ(a)⁻D⁺_yvⁿ_I,j . (3.9) At the boundaries y= 0 andy= 1, it is assumed that

f₀ⁿ(a)≤0, f_Jⁿ(a)≥0, so for j= 0 andj=J, the scheme (3.7) reads

v_i,1ⁿ⁺¹−v_i,0ⁿ

∆τ −σiD⁺_x(D⁻_xv_i,1ⁿ⁺¹)−µ⁺_i D⁺_xvⁿ⁺¹_i,1 −µ⁻_i D⁻_xvⁿ⁺¹_i,1

= max

a∈A(yj) rⁿ_i,1(a)−f_jⁿ(a)D⁺_yv_i,1ⁿ (3.10) and

v_i,Jⁿ⁺¹−vⁿ_i,J

∆τ −σ_iD⁺_x(D⁻_xv_i,Jⁿ⁺¹)−µ⁺_i D⁺_xvⁿ⁺¹_i,J −µ⁻_i D⁻_xv_i,Jⁿ⁺¹

= max

a∈A(yj)

rⁿ_i,J(a)−f_jⁿ(a)D⁻_yv_i,Jⁿ (3.11)

3.4 Monotonicity

In this section we will show that, provided a linear CFL-condition (3.13), the numerical scheme (3.7) ismonotone. That is, ifuandvare solutions of (3.7), then v⁰≥u⁰ =⇒ vⁿ≥uⁿ ∀n >0, (3.12) provided that

∆τkfk_∞≤∆y. (3.13)

Remark 3.4. The inequalities in (3.12) are to be intended component vice.

(28)

The proof is by induction, suppose thatvⁿ≥uⁿ and choosei^∗, j^∗ such that (i^∗, j^∗) = arg min

i,j (v_i,jⁿ⁺¹−uⁿ⁺¹_i,j ). (3.14) The numerical scheme (3.7) can be rewritten as

vⁿ⁺¹_i,j −∆τ(σ_iD⁺_x(D⁻_xv_i) +µ⁺_i D⁺_xv_i+µ⁻_i D⁻_xv_i) =v_i,jⁿ + ˜H(vⁿ_i,j), (3.15) with

H(v˜ ⁿ_i,j) = max

a∈A(y_j) r_ijⁿ(a)−f_jⁿ(a)⁺D⁻vⁿ_i,j−f_jⁿ(a)⁻D⁺_yv_i,jⁿ .

By writing out all the terms in equation (3.15), we get v_i,jⁿ⁺¹−∆τ σ_iv_i+1,jⁿ⁺¹ −2v_i,jⁿ⁺¹+vⁿ⁺¹_i−1,j

∆x² +µ⁺_i v_i+1,jⁿ⁺¹ −vⁿ⁺¹_i,j

∆x +µ⁻_i vⁿ⁺¹_i,j −v_i−1,jⁿ⁺¹

∆x

!

=v_i,jⁿ + ˜H(v_i,jⁿ ), which can be rewritten as

v_i,jⁿ⁺¹ ∆τ

∆x ∆x

∆τ +2σi

∆x+µ⁺_i −µ⁻_i

−v_i+1,jⁿ⁺¹ ∆τ

∆x σi

∆x+µ⁺_i

−vⁿ⁺¹_i−1,j ∆τ

∆x σ_i

∆x−µ⁻_i

=vⁿ_i,j+ ˜H(v_i,jⁿ ).

(3.16) By rearranging the terms in the last equation, we get

vⁿ⁺¹_i,j =c₁ v_i,jⁿ +c₂v_i+1,jⁿ⁺¹ +c₃vⁿ⁺¹_i−1,j

+c₁H(vⁿ_i,j) (3.17) with

c1=

∆x

∆τ

∆x

∆τ +^2σ_∆xⁱ +µ⁺_i −µ⁻_i , c2= ∆τ

∆x σi

∆x+µ⁺_i

, c3=∆τ

∆x σi

∆x−µ⁻_i , (3.18) and it can be easily checked that c₁, c₂ andc₂ are non-negative. We observe that H(v˜ ⁿ_i,j) can be written as

H(v˜ _i,jⁿ ) = max

rⁿ_ij(a)− 1

∆y

v_i,jⁿ |f_jⁿ(a)|+vⁿ_i,j−1f_jⁿ(a)⁺−v_i,j+1ⁿ f_jⁿ(a)⁻

,

hence

H(v˜ ⁿ_i,j)−H(u˜ ⁿ_i,j)≥min

H(v˜ _i,jⁿ )−H(u˜ ⁿ_i,j)

≥ 1

∆y min

−(v_i,jⁿ −uⁿ_i,j)|f_jⁿ(a)|+ (v_i,jⁿ −uⁿ_i,j−1)f_jⁿ(a)⁺

−(v_i,j+1ⁿ −uⁿ_i,j+1)f_jⁿ(a)⁻

. (3.19)

(29)

From the previous inequality and the fact that vⁿ≥uⁿ, we get H(v˜ ⁿ_i,j)−H(u˜ ⁿ_i,j)≥ − 1

∆y(vⁿ_i,j−uⁿ_i,j)|f_jⁿ|∞, (3.20) with

|f_jⁿ|_∞:= max

a (f_jⁿ(a)).

From equation (3.17) and inequality (3.20), it follows that

v_i,jⁿ⁺¹−uⁿ⁺¹_i,j =c₁ v_i,jⁿ −uⁿ_i,j+c₂(vⁿ⁺¹_i+1,j−uⁿ⁺¹_i+1,j) +c₃(v_i−1,jⁿ⁺¹ −uⁿ⁺¹_i−1,j) +c1∆τ(H(v_i,jⁿ )−H(v_i,jⁿ ))

≥c1 v_i,jⁿ −uⁿ_i,j+c2(vⁿ⁺¹_i+1,j−uⁿ⁺¹_i+1,j) +c3(v_i−1,jⁿ⁺¹ −uⁿ⁺¹_i−1,j)

−c₁∆τ

∆y(vⁿ_i,j−uⁿ_i,j)kf_jⁿk_∞

=c1 c2(v_i+1,jⁿ⁺¹ −uⁿ⁺¹_i+1,j) +c3(vⁿ⁺¹_i−1,j−uⁿ⁺¹_i−1,j) +c1(v_i,jⁿ −uⁿ_i,j)(1−∆τ

∆ykf_jⁿk_∞), so provided that

∆τkf_jⁿk_∞≤∆y, (3.21)

the following inequality holds

v_i,jⁿ⁺¹−uⁿ⁺¹_i,j ≥c₁c₂(v_i+1,jⁿ⁺¹ −uⁿ⁺¹_i+1,j) +c₁c₃(v_i−1,jⁿ⁺¹ −uⁿ⁺¹_i−1,j).

The previous inequality together with (3.14) implies that

v_iⁿ⁺¹∗,j^∗−uⁿ⁺¹_i∗,j^∗ ≥c1c2(vⁿ⁺¹_i∗,j^∗−uⁿ⁺¹_i∗,j^∗) +c1c3(v_iⁿ⁺¹∗,j^∗−uⁿ⁺¹_i∗,j^∗), which can be rearranged into

(1−c1(c2+c3))vⁿ⁺¹_i∗,j^∗ ≥(1−c1(c2+c3))uⁿ⁺¹_i∗,j^∗. From the definition of c1, c2 andc3, given by (3.18), it follows that

(1−c1(c2+c3))>0, this concludes the proof.

3.5 Stability

Let v represent a solution of (3.7). Provided that the CFL-condition (3.13) is satisfied, there existsC >0, independent of ∆x,∆y and ∆τ, such that

kvⁿk∞≤C kv⁰k∞+krk∞

, ∀n >0. (3.22)

(30)

The proof is by induction. Let (i, j) = arg max

i,j (vⁿ_i,j), equation (3.16) implies that

vⁿ⁺¹_i,j ≤vⁿ_i,j+ ∆τH(v˜ ⁿ_i,j).

By writing out the term ˜H(vⁿ_i,j) in the last inequality, we get vⁿ⁺¹_i,j ≤vⁿ_i,j+ ∆τmax

a

rⁿ_ij(a)− 1

∆y

v_i,jⁿ |f_jⁿ(a)|+vⁿ_i,j−1f_jⁿ(a)⁺−v_i,j+1ⁿ f_jⁿ(a)⁻

=vⁿ_i,j+ ∆τ r_ijⁿ(a^∗)−∆τ

∆y

v_i,jⁿ |f_jⁿ(a^∗)|+v_i,j−1ⁿ f_jⁿ(a^∗)⁺−vⁿ_i,j+1f_jⁿ(a^∗)⁻

=vⁿ_i,j(1−∆τ

∆y|f_jⁿ(a^∗)|) + ∆τ vⁿ_i,j−1f_jⁿ(a^∗)⁺−vⁿ_i,j+1f_jⁿ(a^∗)⁻

+ ∆τ r_ijⁿ(a^∗).

(3.23) Suppose that (3.13) is satisfied, since f_jⁿ(a^∗)⁺ ≥0, −f_jⁿ(a^∗)⁻ ≥ 0 and we have choseni, j such thatv_i,jⁿ⁺¹= max

i,j vⁿ⁺¹_i,j , inequality (3.23) implies that max

i,j (vⁿ⁺¹_i,j )≤∆τ

kfⁿk_∞max

i,j (v_i,jⁿ ) +krⁿk_∞

.

By choosing (i, j) = arg min

i,j (v_i,jⁿ ) it follows by similar arguments that min

i,j (vⁿ⁺¹_i,j )≥∆τ

kfⁿk_∞min

i,j (v_i,jⁿ ) +krⁿk_∞

.

The last two inequalities together implies that

kvⁿ⁺¹k_∞≤∆τ(kfⁿk_∞|vⁿ|_∞+krⁿk_∞)

≤∆y|vⁿ|∞+ ∆τkrⁿk∞,

where we have used (3.13) to obtain the last inequality. From the last inequality it follows that

kvⁿ⁺¹k_∞≤∆y(∆ykvⁿ⁻¹k_∞+ ∆τkrⁿ⁻¹k_∞) + ∆τkτⁿk_∞

= ∆y²kvⁿ⁻¹k∞+ ∆y∆τkrⁿ⁻¹k∞+ ∆τkrⁿk∞

≤∆y^k+1kv⁰k∞+

n

X

k=0

∆y^k∆τkr^n−kk, since ∆y≤y_max:= 1, the last inequality implies that

kvⁿ⁺¹k∞≤C kv⁰k∞+krk∞ , withC= max(T,1). This concludes the proof.

(31)

3.6 Consistency

We prove that the scheme given by (3.7) is consistent, provided that the exact solution of (3.6), denotedv(x, y, τ), is sufficiently smooth.

Suppose thatv is second order differentiable inx, first order differentiable iny andτ, then Taylor’s theorem implies that

∂v

∂τ(xi, yj, τn) = v_i,jⁿ −vⁿ⁻¹_i,j

∆τ +O(∆τ), ∂²v

∂x²(xi, yj, τn) = D⁺_x(D⁻_xv_i,jⁿ ) +O(∆x²),

∂v

∂x(xi, yj, τn) = D⁻_xvⁿ_i,j+O(∆x), ∂v

∂x(xi, yj, τn) = D⁺_xv_i,jⁿ +O(∆x),

∂v

∂y(x_i, y_j, τ_n) = D⁺_yvⁿ_i,j+O(∆y), ∂v

∂y(x_i, y_j, τ_n) = D⁻_yv_i,jⁿ +O(∆y).

(3.24) Let the functional H(·) be defined such that

H(v)(x, y, τ) = sup

a∈A(y)

r(x, y, a, τ)−f(y, a, τ)∂v

∂y(x, y, t)

.

Becausev satisfies the equation ∂v

∂τ −1 2σ²∂²v

∂x² −µ∂v

∂x+ρv

(x, y, τ)−H(v)(x, y, τ) = 0,

there holds Z τ_n+1

τn

∂v

∂τ −1 2σ²∂²v

∂x² −µ∂v

∂x +ρv

(x, y, τ)dτ− Z τ_n+1

τn

H(v)(x, y, τ)dτ = 0.

By approximating the two integrals in the last equation with the right end point rule and the left endpoint rule, respectively, we get

v(x, y, τn+1)−v(x, y, τn+1)−∆τ 1

2σ²∂²v

∂x²+µ∂v

∂x −ρv

(x, y, τn+1)

−∆τH(v)(x, y, τn) =O(∆τ²), which can be rewritten as

v(x, y, τn+1)−v(x, y, τn+1)

∆τ −

1 2σ²∂²v

∂x² +µ∂v

∂x−ρv

(x, y, τn+1)

−H(v)(x, y, τn) =O(∆τ).

(32)

From the previous equation and (3.24), it follows that eⁿ⁺¹_i,j :=

v(x_i, y_j, τ_n+1)−v(x_i, y_j, τ_n+1)

∆τ −

1 2σ²∂²v

∂x² +µ∂v

∂x−ρv

(x_i, y_j, τ_n+1)

−H(v)(xi, yj, τn)

−

v_i,jⁿ⁺¹−vⁿ_i,j

∆τ −σiD⁺_x(D⁻_xv_i,iⁿ⁺¹)−µ⁺_i D⁺_xvⁿ⁺¹_i,j

−µ⁻_i D⁻_xvⁿ⁺¹_i,j −H(v˜ ⁿ_i,j)

= ˜H(v_i,jⁿ )−H(v)(xi, yj, τn) +O(∆x+ ∆τ).

We have

H(v˜ ⁿ_i,j)−H(v)(xi, yj, τn) ≤max

a∈Aj

H(v˜ ⁿ_i,j)−H(v)(xi, yj, τn)H

≤max

a∈Aj

−f_jⁿ(a)⁺D⁻v_i,jⁿ −f_jⁿ(a)⁻D⁺vⁿ_i,j +f_jⁿ(a)∂v

∂y(xi, yj, τn)

= max

a∈Aj

|f_jⁿ(a)|O(∆y)

=O(∆y).

This concludes the proof.

3.7 Implementation

The numerical scheme given by (3.7) can be rewritten in matrix form as follows 1

∆τ(vⁿ⁺¹−vⁿ) + Avⁿ⁺¹=rⁿ+ Cⁿvⁿ, n= 0, . . . , N−1, (3.25) where v⁰ is a given initial solution. Before we define the vectorsvⁿ⁺¹,vⁿ and rⁿ and the matrices A and Cⁿ, appearing in the last equation, let

m=I·J, and define the functionξ(x, y, τ), such that

ξ_i,jⁿ = arg sup r_ijⁿ(a)−f_jⁿ(a)⁺D⁻vⁿ_i,j−f_jⁿ(a)⁻D⁺_yv_i,jⁿ

. (3.26)

In equation (3.25) we have introduced the solution vector vⁿ ∈ R^m, for n = 0, . . . , N, given as

vⁿ= (v_1,1ⁿ , . . . v_I,1ⁿ , v_1,2ⁿ , . . . v_I,2ⁿ , . . . , v_1,Jⁿ . . . vⁿ_I,J)^>, the vectorrⁿ∈R^m, forn= 0, . . . , N, given as

rⁿ= (r^∗,n_1,1, . . . r_I,1^∗,n, r^∗,n_1,2, . . . r_I,2^∗,n, . . . , r^∗,n_1,J. . . r^∗,n_I,J)^>, r^∗,n_i,j =r(xi, yj, ξ_i,jⁿ , τn),

Numerical Methods for Valuation and Optimal Operation of Natural Gas Storage