Predicting demand of natural gas: a system identification approach

Predicting demand of natural gas: a system identification approach

Master research project Applied Mathematics

March 2013

Student: R.R.A Prins

First advisor: Prof. Dr. M.K. Camlibel

Abstract

In this thesis we construct a linear discrete-time time-invariant input/state/output sys-

tem in order to predict the demand of natural gas based on the available data. A Python

program is written to create this model based on the known input and output. Further-

more a methodology is established to give the best parameters of the system which lead

to the best predictions.

Acknowledgements

I would like to thank GasTerra B.V. for giving me the opportunity to do my final years project at their company. Furthermore I thank the coworkers at GasTerra B.V. for helping me during the process. I would like to especially than H. Kremer for providing an interesting research, L.M.H. Haveman and S. Douma for their guidance and moral support.

At the university I would like to thank Prof.dr. J. Top for accepting the role as second supervisor. Last but definitely not least I would like to give special thanks to Prof. dr.

M.K. Camlibel for supervising me in this for us both unknown field of research, and for his guidance and support along the way.

Of course I am also grateful for all the support given by friends, family and fellow students, in especially S. Baars for helping me with understanding Python.

Thank you all.

List of Tables

1.1 List of connections to Germany . . . . 2

1.2 List of factors . . . . 3

1.3 Used factors . . . . 3

2.1 Input data . . . . 10

2.2 Output data . . . . 10

3.1 Results for hour-ahead . . . . 20

3.2 Results for 12-hours ahead . . . . 20

3.3 Results for day-ahead . . . . 21

3.4 Results for week-ahead . . . . 21

List of Figures

2.1 n = 17, N = 2500, p = 31, i = 2530 . . . . 18

3.1 Reinitializing . . . . 23

3.2 4,5 years of data . . . . 24

Contents

1 Introduction 1

1.1 GasTerra B.V. . . . . 1

1.2 This research . . . . 1

1.3 Existing methods . . . . 2

1.4 Choosing the model . . . . 3

2 The model 5 2.1 Discrete-time dynamical linear models . . . . 5

2.2 System identification . . . . 7

2.2.1 Preliminaries . . . . 7

2.3 Predicting gas demand . . . . 16

3 Using the model 19 3.1 Reinitialization . . . . 19

3.2 Parameter analysis . . . . 19

4 Conclusions and future work 25 A Python programs 28 A.1 Data processing . . . . 28

A.1.1 Creating the input . . . . 28

A.1.2 Creating the output . . . . 37

A.1.3 Creating input and output . . . . 38

A.2 Estimating the system matrices . . . . 41

A.3 Running the system . . . . 46

A.4 Data analyzis . . . . 47

A.5 Plotting and reinitialization . . . . 49

A.6 Optimizing p,N,n . . . . 54

Chapter 1

Introduction

This master thesis is the result of a research worth 50 European Credits done at GasTerra B.V. in association with the University of Groningen concerns research on the prediction of demand of natural gas.

1.1 GasTerra B.V.

GasTerra B.V. is an international company trading in natural gas. It operates on the European energy market and accounts for a significant share of the Dutch gas supply.

The company also provides services related to gas trading. The company has a strong purchasing position and almost fifty years of experience in natural gas procurement and sales [4].

At GasTerra B.V. there is a need to predict the demand of natural gas, based on fore- casting they can make decisions important to insure their position in the gas market.

The methods used at this moment are based on stationary models, which do not change in response to the alternating environment. Their interest lays at constructing a model which can learn from new information and adjust itself in such a way that it captures the change in response to changing environment.

1.2 This research

Based on the need of GasTerra B.V. we want to construct a model which learns from

the changing environment and predicts the demand of gas. The focus of the thesis will

be on the export to Germany. Table 1.1 shows the connections for pipelines between the

Netherlands and Germany. The connection Zevenaar/Winterswijk will be considered in

this work. This connection is together with Oude StatenZijl (OSZ) the biggest connection

to Germany. Zevenaar/Winterswijk delivers only L-gas, OSZ on the other hand delivers

both L-gas and G-gas. Tegelen also delivers to to the same distributer in Germany, but

in comparison to Zevenaar/Winterswijk this is a really small quantity, therefore it is

not considered in this thesis. Hence we focus on the total amount of natural gas which

Oude Statenzijl (OSZ) Winterswijk

Zevenaar Tegelen Bocholtz

Table 1.1: List of connections to Germany

flows through Zevenaar/Winterswijk. The aim of this work is to create a dynamic model which can predict the demand of natural gas and adjusts itself in case of changes in the behavior of the system.

The rest of the chapter will proceed as follows. First we will give a summary of the already existing methods in Section 1.3. In Section 1.4 the proposed model will be introduced.

1.3 Existing methods

In [3],the authors proposed a model for the demand of gas for household use in four states of the United States of America, first by using only physical factors and second by also implementing economical factors. It turns out that for the four states the model is differently accurate. Furthermore, the influence of economical factors was in some cases insignificant and in other cases significant but not of big influence in comparison to the physical factors. One can conclude from this that physical factors are the leading influence on the demand of natural gas and economical factors might be of influence (but of smaller influence then physical factors). H.G. Berrisford [2] stated that using predetermined seasonal pattern would cause less power to predict. He only uses the direct past history of the temperature by means of an exponentially weighted average.

In [1], P. Balestra and M. Nerlove propose a dynamical model, for the estimation of the parameters they use the least squares and maximum likelihood method, both methods of estimation result in weak predictions. In [17], Siemek et al. employ a Hubert model which is based on a decay method. This method works too slow to be able to deal with hourly data. It is used to estimate an overall demand of natural gas of a specific country, so not specifically for a company, or for a certain client. Liu and Lin [12] forecast household use of natural gas by using monthly and quarterly time series. They found that using quarterly data was less time-consuming, but not necessarily more accurate in comparison to monthly data. Furthermore, they state that using monthly data has the tendency to over-parametrize the system. In [18], the authors model the demand of gas in Bangladesh. Their model gives a lower bound and does not estimate the actual amount.

Many other models exist for predicting the demand of gas. Most of these models are

stationary. Almost all papers estimate the demand of household use, and are not aiming

to estimate the amount of gas for a specific company. In this thesis, we will propose a

dynamical model in order to predict the demand of natural gas for a specific company.

1.4 Choosing the model

In [13], Ljung mentions that the construction of a model involves three basic entities:

1. A data set.

2. A set of candidate models

3. A rule by which the candidate models can be assessed using the data.

This can be used as a guide to decide which kind of model one should use. We start with constructing a data set. The list of factors in Table 1.2 consists of factors which might be of interest. Among these possible factors, we focus n physical parameters suggested by the experts of GasTerra B.V.. The output will be the physical flow.

Temperature Price of gas Day

Hours sun Wind Price other energy carriers

CO

Price Spark spread Storage

Unemployment rate Consumer trust Einkaufsmanagers index

Transportation costs rainfall Time

Table 1.2: List of factors

Spark spread shows a balance between the price of gas and the price of coal included the CO

price. CO

Price represents the costs a company has to pay for the CO

created by for instance burning coal. Einkaufsmanagers index is an index which stands for the development of the economy of Germany.

For data concerning weather we will retrieve our information from the Koninklijk Ned- erlands Meteorologisch Instituut (KNMI) [9]. For the data concerning gas flows we will retrieve information from the Gasunie Transport Services (GTS) [16]. In the rest of the thesis we will use the factors given in Table 1.3.

Direction of the wind Speed of the wind Temperature Sunshine duration Solar radiation Duration of rainfall Amount of rainfall Air pressure Overcast

Table 1.3: Used factors

The direction of the wind (in degrees) represents the average of the last 10 minutes of the past hour. Speed of the wind (in m/s) is the hourly average. Temperature (in degrees Celsius) is measured 1.50 meter above sea level and represents the temperature at that moment. Sunshine duration (in hours) is the duration of sunshine in an hour.

Solar radiation (in J/cm

) represents the strength of the sun. The duration of rainfall (in hours) is the duration of rainfall during the hour. The amount of rainfall (in mm) is the total amount of rainfall in an hour. Air pressure is given in hPa. Overcast is a number which represents the cloudiness.

In this thesis, we will construct a discrete-time dynamical model based on the factors

given in Table 1.3. One of the advantages working with the model class we have chosen is the case of including new factors and/or discarding existing ones.

The class of discrete dynamical linear models can be divided into time-invariant and time-varying models. The demand of gas is most likely time-varying, since the same conditions in the winter might lead to another result than the same conditions in the summer. However, dealing with time-varying models is not always preferable due to the heavy computations required by such models. In this thesis, we will employ the-invariant models and try to capture the time dependency of the model by creating new models over time. The findings of this thesis can be taken as a starting point for further research on modeling the gas demand by time-varying systems. We will discuss the details of the model we propose in Section 2.1. Before dealing with the mathematical details, we first discuss various existing methodologies in what follows.

In [6], the authors investigate how a linear input/state/output system can be constructed from given input/output data. Their method uses the so-called information matrix along with its Toeplitz matrix structure. The information matrix contains the given input/out- put data. By using the information matrix as well as its shifted versions, they construct the matrices that determine a linear input/state/output system. In [19], Wang and Qin proposed a subspace identification method that is based on principal component analy- sis (PCA) which gives consistent estimates in the error-in-variables (EIV) formulation.

Lacy and Bernstein presented a constraint optimization method to ensure asymptotic stability of the identified model in the context of subspace identification in [11]. Several other papers investigated how to estimate the system matrices in the presence of noise, e.g. [5] and [15]. Katayama and Tanaka [8] also use the noise term, they estimate the matrices using an LQ-decomposition which can be used to create the orthogonal and oblique projections of parts of the information matrices. The paper [10] exploits the information and its Toeplitz structure in the presence of noise. In [7] adaptive Iterative learning control (ILC) has been employed. They first estimate the matrices, in what- ever way one might prefer, and then they keep reshaping the system matrices by use of ILC. Merc` ere and Bako [14] revise the problem of determining a canonical state-space representation for multi-variable systems. They first estimate the system matrices and finally put them in canonical representation.

Most of the methods mentioned in the papers are using an information matrix and do basically the same as Juang did in [6], even the papers concerning noise terms often use this method of estimation. In all the papers mentioned above, the estimation problem boils down to solving certain linear equations. These papers differ from each other by the methods they proposed in solving the equations. Such as LQ-factorization, least squares, or pseudo-inverse etc. In this thesis, we will focus on the method proposed by Juang in [6].

The structure of the thesis is as follows. In Chapter 2 the modeling methodology we

adopt will be discussed in detail. Chapter 3 is devoted to the question of how to update

the model when new data is available. Finally, the thesis closes with the conclusions and

a discussion for further research in Chapter 4. The computer program that was used for

this thesis is included in Appendix A.

Chapter 2

The model

Physical behavior can be modeled in a system, for instance the motion of a car driving in a straight line. We can use as input the pressure given on the throttle, as output we want the distance traveled. This is the best modeled by a dynamical input/state/output system, where the state represents some physical quantities and is able to keep track of the speed and the distance. This chapter will provide a dynamical system which will be used to predict the amount of gas. Note that we do not know the exact behavior of the demand, and therefore cannot give physical meaning to the state.

The organization of this chapter is as follows. In Section 2.1 we will introduce discrete time-invariant input/state/output dynamical linear systems and review certain proper- ties of such systems. In Section 2.2, we will present a methodology based on [6] in order to estimate system matrices from given input/output observations. Finally we will dis- cuss how this methodology can be employed in the context of predicting gas demand in Section 2.3.

2.1 Discrete-time dynamical linear models

Consider the linear discrete-time time-invariant input/state/output system Σ(A, B, C, D) of the form:

x(k + 1) = Ax(k) + Bu(k),

y(k) = Cx(k) + Du(k), (2.1)

where k ∈ N denotes the time, u(k) ∈ R

the input vector, x(k) ∈ R

the state vec-

tor, y(k) ∈ R

the output vector and the so-called system matrices A, B, C, D are of

appropriate dimensions. We sometimes refer to this representation as state-space repre-

sentation. If the system description (2.1) is obtained from a physical system the state

of the system (very often) has a physical interpretation. For instance, if one models an

electrical circuit containing linear dynamical elements like capacitors and inductors then

the components of the state vector correspond to the voltages/charges of the capacitors

and the current/flows of the inductors. In this thesis we attempt to identify a system

of the form (2.1) in order the model gas demand. As such, the state vector will not

necessarily admit a certain physical meaning.

Next, we illustrate the underlying dynamics of the system of the form (2.1) by means of an example.

Example Consider the state space representation of the form (2.1) where A = 1 2

3 4



, B = 1 −0.2 3 4 −0.5 6



, C = 1 2 , D = 1 −0.4 2 .

For this model, the input vector u consists of three components: the first corresponds to the speed of the wind (10

m/s), the second to the temperature (10

‰), and the third corresponds to the amount of rain fallen in the hour (10

mm)

whereas the scalar output y corresponds to the amount of gas (10

kWh)

.

Suppose that the initial state x(0) is given by x(0) = 1

1

 .

To calculate the output we use the input, the system matrices and the state together with equation (2.1). Assume that at time instance 0 the speed of the wind is 1.0 m/s, the temperature is 10.6 ‰ and the amount of rain fallen is 1.6 mm. This leads to the output at time instance 0

y(0) = Cx(0) + Du(0) = 1 2 1 1



+ 1 −0.4 2



 10 106

16



 = 2.6.

Additionally, we calculate x(1):

x(1) = Ax(0) + Bu(0) = 1 2 3 4

 1 1



+ 1 −0.2 3 4 −0.5 6





 10 106

16



 = 39, 8 96

 .

Assume at time instance 1 that the speed of the wind is 2.0 m/s, the temperature is 10.8

‰ and the amount of rain fallen is 1.2 mm. This leads to the output at time instance 1

y(1) = Cx(1) + Du(1) = 1 2 39, 8 96



+ 1 −0.4 2



 20 108

12



 = 232.6.

Additionally, we calculate x(2):

x(2) = Ax(1) + Bu(1) = 1 2 3 4

 39, 8 96



+ 1 −0.2 3 4 −0.5 6





 20 108

12



 = 266.2 601.4

 .

In a similar fashion, one can compute the state and the output at time instant k, i.e.

x(k) and y(k), based on the values of inputs u(0), u(1), . . . , u(k).

11 millimeter rainfall is the equivalent of 1 liter of water per square meter per hour.

2.2 System identification

Roughly speaking, the system identification problem we deal with in this thesis amounts to construct a dynamical system of the form (2.1) based on given input/output observa- tions. In order to elaborate on the method we use, we first introduce a set of machinery in what follows

2.2.1 Preliminaries

This section is devoted to the definitions of mathematical concepts and notations that will be used throughout the thesis.

Most of the time we follow the standard mathematical notation and conventions. We denote the set of natural numbers (including zero) by N, the set of real number by R, the set of q-tuples of real numbers by R

and the set of q × p matrices by R

.

In this Chapter the Moore-Penrose pseudo-inverse will be used, hence we define it below.

Definition 2.2.1 (Moore-Penrose pseudo-inverse). For the matrix A ∈ R

, the Moore-Penrose pseudo-inverse is defined as A

∈ R

and satisfies the following prop- erties:

(i) AA

A = A (ii) A

AA

= A

(iii) (AA

)

= AA

(iv) (A

A)

= A

A.

Furthermore, A

exists and is unique.

In the case that A is invertible this inverse will be the normal inverse of an invertible matrix. The advantage of using this pseudo-inverse is that it is widely used and most mathematical programming languages have a build in program to compute the Moore- Penrose pseudo-inverse of any matrix.

A commonly used class of matrices in the context of system identification is the so-called Hankel matrices. The Output matrix and the Input matrix have the shape of a Hankel matrix.

Definition 2.2.2 (Hankel matrix). A matrix M ∈ R

is called a Hankel matrix if M

= M

whenever i + j = k + l, that is if it is off the form:















m

m

m

· · · m

m

m

m

· · · m

m

m

m

· · · m

.. . .. . .. . . .. .. . m

m

m

· · · m













 .

for some real numbers m

with 1 ≤ i ≤ 2q − 1.

We sometimes say that a (possibly non-square) matrix is a block Hankel matrix if it admits the Hankel structure blockwise.

Another relevant class of matrices is the so called Toeplitz matrices.

Definition 2.2.3 (Toeplitz matrix). Given the system matrices A ∈ R

, B ∈ R

, C ∈ R

, D ∈ R

and an integer p ∈ N we define the corresponding Toeplitz matrix as

T

:=





















D 0

CB D 0

CAB CB D 0

.. . .. . .. . . .. ...

.. . .. . .. . .. . . .. 0 CA

B CA

B CA

B · · · · · · D



















 .

In the course of the discussion, we will work with matrices obtained from sequences of inputs, states and outputs. In what follows we introduce these matrices.

Definition 2.2.4 (Input matrix). Let u(k) ∈ R

be the input at time instance k ∈ N and p ∈ N, then we define the vector u

(k) ∈ R

as:

u

(k) :=















u(k) u(k + 1) u(k + 2)

.. . u(k + p − 1)













 .

Then for given N ∈ N we define the input matrix U

(k) ∈ R

as:

U

(k) := u

(k) u

(k + 1) · · · u

(k + N − 1) 

=











u(k) u(k + 1) · · · u(k + N − 1) u(k + 1) u(k + 2) · · · u(k + N )

.. . .. . . .. .. .

u(k + p − 1) u(k + p) · · · u(k + p + N − 2)









 .

Similar to the input matrix, we define the output matrix as the following.

Definition 2.2.5 (Output matrix). Let y(k) ∈ R

be the output at time instance k ∈ N and p ∈ N, then we define the vector y

(k) ∈ R

as:

y

(k) :=















y(k) y(k + 1) y(k + 2)

.. .















.

Then for given N ∈ N we define the output matrix Y

(k) ∈ R

as:

Y

(k) := y

(k) y

(k + 1) · · · y

(k + N − 1) 

=











y(k) y(k + 1) · · · y(k + N − 1) y(k + 1) y(k + 2) · · · y(k + N )

.. . .. . . .. .. .

y(k + p − 1) y(k + p) · · · y(k + p + N − 2)









 .

Note that the input and output matrix are of block Hankel matrix form.

Apart from these matrices, we need to define the state matrix.

Definition 2.2.6 (State matrix). Let x(k) ∈ R

be the state at time instance k ∈ N then for a given N ∈ N we define the state matrix X(k) ∈ R

as:

X(k) := x(k) x(k + 1) · · · x(k + N − 1) . Finally, we define the observability matrix.

Definition 2.2.7 (Observability matrix). Given the system matrices A ∈ R

, C ∈ R

and p ∈ N we define the observability matrix as:

O

:=



















 C CA CA

.. . .. . CA



















 .

The relationships among the matrices defined above is in order. By straightforward calculations, we arrive at the following state and output trajectories of the system (2.1) for the initial state x(0) and the input sequence u(0), u(1), . . . , u(k):

x(k) = A

x(0) + A

Bu(0) + A

Bu(1) + · · · + ABu(k − 2) + Bu(k − 1), y(k) = CA

x(0) + CA

Bu(0) + · · · + CABu(k − 2) + CBu(k − 1) + Du(k).

In a more compact form we have

y

(k) = O

x(k) + T

u

(k), (2.2)

Y

(k) = O

X(k) + T

U

(k). (2.3)

The goal of the method we use it to determine A, B, C, D matrices from given input/out-

put observations. To do so, we first focus on determining the observabilitymatrix O

.

Define:

R

= 1

N Y

(k)Y

(k), R

= 1

N Y

(k)U

(k), R

= 1

N U

(k)U

(k), R

= 1

N X(k)X

(k), R

= 1

N Y

(k)X

(k), R

= 1

N X(k)U

(k).



 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(2.4)

The matrices R

∈ R

, R

∈ R

, R

∈ R

are symmetric matrices that represent the autocorrelations of the output data y with the time shifts, the input data u with the time shifts, and the state vector x, respectively, whereas the matrices R

∈ R

, R

∈ R

, R

∈ R

represent the cross correlations of the output data y with time shifts and the input data u with time shifts, the output data y with time shifts and state vector x, and the state vector x and the input u with time shifts, respectively.

0 1 2 3 4 5

Speed of the wind (in 10

m/s) 10 10 20 30 10 20 Temperature (in 10

‰) 98 94 136 106 102 92 Amount of rain fallen (in 10

mm) 18 73 16 33 40 31

Table 2.1: Input data

0 1 2 3 4 5

Energy (in 10

kWh) 1.2 1.1 1.2 1.1 1.2 1.3 Table 2.2: Output data

Example We can use the input data and output data given in Table 2.1 and Table 2.2,

respectively to calculate the input matrix U

(0) and the output matrix Y

(0). We have

r = 3 and m = 1, choose p = 3 and N = 4 and k = 0. This leads to the following input

and output matrices:

U

(0) =

























10 10 20 30 98 94 136 106 18 73 16 33 10 20 30 10 94 136 106 102 73 16 33 40 20 30 10 20 136 106 102 92 16 33 40 31























 ,

Y

(0) =

" 1.2 1.1 1.2 1.1 1.1 1.2 1.1 1.2 1.2 1.1 1.2 1.3

# .

Using these matrices the quantities R

, R

, R

can be calculated using (2.4) Estimating the observability matrix

By post-multiplying (2.3) by U

(k)/N , we get

R

= Y

(k)U

(k) = O

X(k)U

(k) + T

U

(k)U

(k) = O

R

+ T

R

.

Suppose that R

is nonsingular. Then, we obtain the Toeplitz matrix T

by post- multiplying the last equation by R

:

T

= (R

− O

R

) R

(2.5) Similarly, we get

R

= Y

(k)Y

(k) = O

X(k)Y

(k) + T

U

(k)Y

(k) = O

R

+ T

R

, (2.6) and

R

= Y

(k)X

(k) = O

X(k)X

(k) + T

U

(k)X

(k) = O

R

+ T

R

. (2.7) by post-multiplying (2.3) by Y

/N and by X

/N , respectively. Substituting (2.7) in (2.6) results in

R

= O

O

R

+ T

R

+ T

R

= O

R

O

+ O

R

T

+ T

R

. Furthermore, one can substitute T

given by (2.5) and get:

R

= O

R

O

+ O

R

T

+ T

R

,

= O

R

O

+ O

R

(R

− O

R

) R

+ (R

− O

R

) R

R

,

= O

R

O

+ O

R

R

R

− O

R

R

R

O

+ R

R

R

− O

R

R

R

,

= O

R

O

− O

R

R

R

O

+ R

R

R

.

This leads to

R

− R

R

R

= O

R

− R

R

R

O

. (2.8) Now, we define

R

= R

− R

R

R

, (2.9) R ˜

= R

− R

R

R

. (2.10) Then, (2.8) can be written as

R

= O

R ˜

O

. (2.11)

Clearly, both R

∈ R

and ˜ R

∈ R

are symmetric matrices.

We are now in a position to compute an estimate of the observability matrix. By following the footsteps of [6], let R

R

= U Σ

V

= U

U



 Σ

0

0

0

 V

V



= U

Σ

V

, (2.12) be the singular value decomposition of R

where n

is the amount of singular values equal to 0. Since R

is symmetric, we have U

= V

and U

= V

. Then, it follows from (2.11) that

R

= O

R ˜

O

= U

Σ

U

This leads to the obvious estimates for O

and ˜ R

:

O

= U

, (2.13)

R ˜

= Σ

. (2.14)

In this method the dimension of the state is equal to n. In order to get a specific dimension of the state, we must have pm − n singular values equal to 0. In practice non of the singular values of R

are equal to 0. To achieve the desired dimension of the state we will truncate the smallest pm − n singular values. So Σ

contains the n highest singular values. Then, we have U

∈ R

= O

.

Example Assume p = 3 and N is chosen such that R

is non-singular. Using R

, R

, R

we can calculate R

. Assume that we have

R

=





2 3 5 3 8 7 5 7 4



 .

Using the singular value decomposition shown in (2.12) gives us the following matrices for U

, Σ

, V

:

U

=





−0.39 −0.57 −0.73

−0.71 −0.33 0.63

−0.59 0.75 −0.28



 ,

Σ

=





15.49 0 0

0 2.85 0

0 0 1.36



 ,

V

=





−0.39 −0.57 −0.73

−0.71 −0.33 0.63

−0.59 0.75 −0.28



 . Clearly, U

= V

. Assume n = 2, then this leads to:

O

=





−0.39 −0.57

−0.71 −0.33

−0.59 0.75



 . (2.15)

Estimating A

We construct the matrix A ∈ R

using the observability matrix. We have the follow- ing:

O

(m + 1 : pm, :) =













 CA CA

CA

.. . CA















=













 C CA CA

.. . CA















A = O

(1 : (p − 1)m, :)A.

We use the Moore-Penrose pseudo-inverse to compute the matrix A. This leads to:

A = O

(1 : (p − 1)m, :)O

(m + 1 : pm, :), (2.16) where O

(1 : (p−1)m, :) denotes the Moore-Penrose pseudo-inverse of O

(1 : (p−1)m, :).

Example Using the observability matrix given by (2.15) we have A = −0.39 −0.57

−0.71 −0.33



−0.71 −0.33

−0.59 0.75



= 0.39 −1, 90 0.98 1.83



. (2.17)

Estimating C

The estimation of the matrix C ∈ R

follows directly from the estimation of the observability matrix. Recall that the first m rows of the observability matrix are exactly C. In other words, we have

C = O

(1 : m, :). (2.18)

Example Using the observability matrix given by (2.15) we have

C = −0.39 −0.57 . (2.19)

Estimating B and D

The observability matrix and the form of the Toeplitz matrix are used for the estimation of the matrices B and D. To do so we use the output-error minimization method proposed in [6].

We begin with substituting p = N and k = 0 in (2.2)

y

(0) = O

x(0) + T

u

(0). (2.20) We write the matrices B and D in the following way:

B = b

b

· · · b

 , D = d

d

· · · d

 ,

where b

∈ R

and d

∈ R

for i = 1, 2, . . . , r. Define the vectors b ∈ R

and d ∈ R

as:

b =









 b

b

.. . b











, d =









 d

d

.. . d









 .

Also, let

u(k) =









 u

(k) u

(k)

.. . u

(k)









 .

We now have the following for T

u

(0):

T

u

(0) =













 D

CB D

CAB CB D

.. . .. . .. . . ..

CA

B CA

B CA

B · · · D



























 u(0) u(1) u(2) .. . u(N − 1)













 ,

=















U

(0) U

(1) U

(2)

.. . U

(N − 1)













 d +



















0

CU

(0) CAU

(0) + CU

(1)

.. .

N −2

P

k=0

CA

U

(k)

















 b,

(2.21)

where 0

∈ R

denotes the zero matrix, and

U

(k) = I

u

(k) I

u

(k) · · · I

u

(k) , U

(k) = I

u

(k) I

u

(k) · · · I

u

(k) ,

where I

∈ R

and I

∈ R

are the identity matrices of order m and n, respec- tively. Substituting the expression for T

u

(0) from (2.21) into (2.20) gives:

y

(0) = ΦΘ, where

Θ =



 x(0)

d b



 , Φ =



















C U

(0) 0

CA U

(1) CU

(0)

CA

U

(2) CAU

(0) + CU

(1)

.. . .. . .. .

CA

U

(N − 1)

N −2

P

k=0

CA

U

(k)

















 .

Using the Moore-Penrose pseudo-inverse, we get the following:

Θ = Φ

y

(0), (2.22)

where Φ

denotes the Moore-Penrose pseudo-inverse of Φ. Then we have:

b = Θ(mr + n + 1 : (r + 1)n + mr), d = Θ(n + 1 : mr + n). (2.23) This method also provides an estimation of the initial state vector x(0):

x(0) = Θ(1 : n). (2.24)

Example Assume p = 1, N = 4, m = 1, r = 3 and A and C are given by (2.17) and (2.19), respectively. Furthermore, assume that we have the input and output data given by Table 2.1 and Table 2.2, respectively. Then, we get

Φ =







−0.39 −0.57 10 98 18 0 0 0 0 0 0

−0.41 −0.30 10 94 73 −3.9 −5.7 −38.22 −55.86 −7.02 −10.26

−0.14 0.22 20 136 16 −7.97 −8.72 −76.50 −83.19 −35.79 −47.05 0.27 0.66 30 106 33 −13.24 −12.23 −104.73 −84.41 −39.39 −27.22





 , and

y

(0) =







 1.2 1.1 1.2 1.1









.

As such, we have

Φ

=







































1.66 · 10

−2.26 · 10

−2.13 · 10

2.92 · 10

−1.46 · 10

−1.40 · 10

−6.46 · 10

1.74 · 10

9.54 · 10

−8.08 · 10

−9.13 · 10

1.12 · 10

1.15 · 10

−2.70 · 10

−1.38 · 10

1.18 · 10

−1.23 · 10

1.92 · 10

5.82 · 10

−1.27 · 10

−2.12 · 10

2.44 · 10

2.75 · 10

−3.78 · 10

5.50 · 10

1.01 · 10

8.51 · 10

−6.88 · 10

5.24 · 10

2.81 · 10

−2.03 · 10

−4.81 · 10

−2.90 · 10

1.02 · 10

1.18 · 10

−1.28 · 10

8.32 · 10

−1.22 · 10

−9.65 · 10

5.57 · 10

1.80 · 10

−9.66 · 10

−1.20 · 10

−1.73 · 10





































 .

Together with (2.22), this leads to

Θ =









































−4.49 · 10

4.09 · 10

5.04 · 10

1.21 · 10

−1.88 · 10

−1.10 · 10

4.73 · 10

1.12 · 10

−3.59 · 10

3.75 · 10

8.01 · 10









































. (2.25)

Using Θ, (2.23) results in

B = −1.10 · 10

1.12 · 10

3.75 · 10

4.73 · 10

−3.59 · 10

8.01 · 10



, (2.26)

D = 5.04 · 10

1.21 · 10

−1.88 · 10

 , (2.27) and (2.24) in

x(0) = −4.49 · 10

4.09 · 10



. (2.28)

2.3 Predicting gas demand

The methodology described in this chapter can be used to construct the linear discrete-

time time-invariant input/state/output system Σ(A, B, C, D) given by (2.1). The ma-

trices A, B, C, D and the initial state are given by (2.16),(2.23),(2.18),(2.23) and (2.24),

Note that the initial state is calculated for the first input/output used for the estima- tions, so u(0) is the first input. For the creation of the matrices we used N + p − 1 data-points. So the output y(N + p − 1), resulting from input u(N + p − 1) will be the first prediction.

Python is used to create programming that can calculate the system matrices A, B, C, D and the initial state x(0), and create a system which can produce the estimated output y for given input u. Additionally, we created programming to plot the actual output against the estimated output as well as the error. The code can be found in Appendix A. In Figure 2.3, such a plot is given. For the creation of the matrices we have used n = 17, N = 2500, p = 30. The parameter i = 2530 indicates the time-index after which only predicted output is estimated, in the plot this is indicated by the vertical line. The second part of the plot contains the estimation error. The horizontal lines indicate a certain allowed error set for this example at 0.25 · 10

. One can see that the error for the first 50 hours of predictions lay within the allowed error range. Later the error grows.

This result implies that the model needs to be reinitialized after a certain amount of

hours to decrease the error. Another point of interest is, how to choose n, N, p. These

points will be further addressed in Chapter 3.

0 1000 2000 3000 4000 5000 6000 7000 8000 Time (Hours)

1 0 1 2 3 4 5

Energy (10^7 kWh)

Estimation Real

0 1000 2000 3000 4000 5000 6000 7000 8000 Time (Hours)

3 2 1 0 1 2 3

Error (10^7 kWh)

Figure 2.1: n = 17, N = 2500, p = 31, i = 2530

Chapter 3

Using the model

In Chapter 2 we discussed the tools to create the system matrices. This chapter aims at developing a method to decide which n, p, N should be chosen in order to obtain a model that provides the best predictions. The structure of this chapter as follows. In Section 3.1 we will discuss how to reinitialize the model. Section 3.2 investigates the system parameters for different length of predictions.

3.1 Reinitialization

The world of natural gas consumption is changing over time, therefore a time-varying model will become less accurate over time. This leads to the idea to reinitialize the system, in such a way that these changes will be captured in the model. We use the error of the system and/or a certain time-step as a criterion to reinitialize. At the moment that the predictions obtained by the model becomes less accurate we would like to reinitialize the model. For the best predictions we would like to use the newest data, hence we would like to reinitialize the system after a certain period of time. Figure 3.1 displays 4 plots from the reinitialization procedure. From i = 2530 on we predict, every subplot shows the new prediction from the given i-value on. On can see that this procedure create better predictions using the reinitialization. Appendix A.5 contains the written code for the reinitialization.

3.2 Parameter analysis

This section provides a method to come up with the best n, N, p for different prediction lengths. We distinguish between 1 hour ahead, 12 hours-ahead, 1 day-ahead, 1 week- ahead and longer ahead.

Hour-ahead

For hour-ahead predictions, we use the following strategy:

n 10 11 12 13 14 15

Best p 19 16 16 15 19 20

Best N 840 840 720 840 720 720

Average quality (2-norm) 0.052 0.046 0.040 0.046 0.061 0.062 Table 3.1: Results for hour-ahead

1. For n = 10, 11, . . . , 15, N = 240, 360, . . . , 1440 and p = n + 2, n + 3, . . . , n + 10 2. Create hour-ahead predictions of the last 30 hours.

3. Calculate the 2-norm of the vector of errors for the predicted hour.

4. Take the average of these 2-norms for given p.

5. Determine the values of n, p and N that give the best results.

For hour-ahead, it turns out that choosing N greater than 840 does not produce better predictions. The best predictions occur when n = 12. This value for n produces the best values for all good p and N . For n = 12 the best value for p fluctuates between 15 and 19. The average prediction error is mostly the same for p between 14 and 17. Whereas the best results are obtained for the values p = 15 and p = 16.

So, for hour-ahead this method suggests to use the following parameters for prediction:

N = 840, n = 12, p = 15 or 16.

12 hours-ahead

n 10 11 12 13 14 15

Best p 13 13 14 15 19 20

Best N 1680 1680 1680 1680 1680 1680 Average quality (2-norm) 0.197 0.194 0.187 0.189 0.222 0.224

Table 3.2: Results for 12-hours ahead For 12 hours-ahead predictions, we use the following strategy:

1. For n = 10, 11, . . . , 15, N = 1200, 1320, . . . , 2400 and p = n + 2, n + 3, . . . , n + 10 2. Create 12 hours-ahead predictions of the last 30 times 12 hours.

3. Calculate the 2-norm of the vector of errors for the predicted 12 hours.

4. Take the average of these 2-norms for given p.

5. Determine the values of n, p and N that give the best results.

For 12 hours-ahead, it turns out that choosing N greater than 1680 does not produce better predictions. The best predictions occur when n = 12, but n = 13 is mostly almost as good. This value for n produces the best values for all good p and N . For n = 12 and n = 13 the best value for p fluctuates between p = 13 and p = 15. Whereas the best results are obtained for the values p = 14 and p = 15.

So for 12 hours-ahead this method suggests to use the following parameters for prediction:

N = 840, n = 12 or 13, p = 14 or 15.

Day-ahead

n 10 11 12 13 14 15

Best p 14 15 17 18 19 20

Best N 1920 2160 1920 1920 1920 2160 Average quality (2-norm) 0.244 0.247 0.222 0.223 0.231 0.239

Table 3.3: Results for day-ahead For day-ahead predictions, we use the following strategy:

1. For n = 10, 11, . . . , 15, N = 1440, 1680, . . . , 3840 and p = n + 2, n + 3, . . . , n + 10 2. Create day-ahead predictions of the last 30 days.

3. Calculate the 2-norm of the vector of errors for the predicted day.

4. Take the average of these 2-norms for given p.

5. Determine the values of n, p and N that give the best results.

For day-ahead, it turns out that choosing N greater than 1920 does not lead to better predictions. On the average n = 12 and n = 13 result in the best predictions. For n = 12 and n = 13 the best value for p changes between 17 and 18. Both values result in practically the same error.

So for day-ahead this method suggests to use the following parameters for prediction:

N = 1920, n = 12 or 13, p = 17 or 18.

Week-ahead

n 10 11 12 13 14 15

Best p 15 17 18 18 19 20

Best N 2736 2736 2736 2736 2736 2736 Average quality (2-norm) 0.346 0.273 0.231 0.217 0.214 0.229

Table 3.4: Results for week-ahead

For week-ahead predictions, we use the following strategy:

1. For n = 10, 11, . . . , 15, N = 1680, 1848, . . . , 3360 and p = n + 2, n + 3, . . . , n + 10 2. Create week-ahead predictions of the last 30 weeks.

3. Calculate the 2-norm of the vector of errors for the predicted week.

4. Take the average of these 2-norms for given p.

5. Determine the values of n, p and N that give the best results.

For week-ahead, it turns out that choosing N greater than 2736 does not lead to better predictions. The best predictions occur when n = 13 and n = 14. Whereas n = 14 predicts slightly better. For n = 13 and n = 14 the best values for p fluctuates between p = 17 and p = 20. Whereas the best results are obtained for the values p = 18 and p = 19.

So for week-ahead this method suggests to use the following parameters for prediction:

N = 2736, n = 13 or 14, p = 17 or 18.

Longer prediction horizons

In Figure 3.2(a) we use one year of data to predict the next 3.5 years, in 3.2(b) 2 years

of data is used to predict the next 2.5 years. So using 2 years of data for predicting

more then a month-ahead leads to better results then 1 year. It is most likely that

using even more data for the estimations result in even better results. But due to lack

of data we are not able to use more data. The input of the system consists of only

weather data, more then a week-ahead prediction of weather data are highly inaccurate

and therefor the predictions given by the system will also be inaccurate. The method

suggested above for the prediction of smaller prediction horizons can be used for longer

prediction horizons, although it does take more time and uses more data. We did take a

look at month ahead, it turns out that n = 12 does lead to the best results together with

p = 17 in the most cases, but finding the correct value for N is too time-consuming.