Performance comparison of the Extended Kalman Filter and the Recursive Prediction Error Method

(1)

Performance comparison of the

Extended Kalman Filter and the

Recursive Prediction Error Method

Examensarbete utfört i Reglerteknik vid Tekniska Högskolan i Linköping

av

Jonas Wiklander Reg nr: LiTH-ISY-EX-3351

(2)

(3)

Performance comparison of the

Extended Kalman Filter and the

Recursive Prediction Error Method

Examensarbete utfört i Reglerteknik vid Tekniska Högskolan i Linköping

av

Jonas Wiklander Reg nr: LiTH-ISY-EX-3351

Supervisors: Geir Hovland (ABB) Mats Larsson (ABB) Markus Gerdin (ISY) Examiner: Svante Gunnarsson Link¨oping 6th June 2003.

(4)

(5)

Avdelning, Institution Division, Department

Institutionen för Systemteknik

581 83 LINKÖPING

Datum Date 2003-06-06 Språk Language Rapporttyp Report category ISBN Svenska/Swedish X Engelska/English Licentiatavhandling

X Examensarbete ISRN LITH-ISY-EX-3351-2003

C-uppsats

D-uppsats Serietitel och serienummer _{Title of series, numbering} ISSN

Övrig rapport

____

URL för elektronisk version

http://www.ep.liu.se/exjobb/isy/2003/3351/ Titel

Title

Jämförelse mellan Extended Kalmanfiltret och den Rekursiva Prediktionsfelsmetoden Performance comparison of the Extended Kalman Filter and the Recursive Prediction Error Method Författare Author Jonas Wiklander Sammanfattning Abstract

In several projects within ABB there is a need of state and parameter estimation for nonlinear dynamic systems. One example is a project investigating optimisation of gas turbine operation. In a gas turbine there are several parameters and states which are not measured, but are crucial for the performance. Such parameters are polytropic efficiencies in compressor and turbine stages, cooling mass flows, friction coefficients and temperatures. Different methods are being tested to solve this problem of system identification or parameter estimation. This thesis describes the implementation of such a method and compares it with previously implemented identification methods. The

comparison is carried out in the context of parameter estimation in gas turbine models, a dynamic load model used in power systems as well as models of other dynamic systems. Both simulated and real plant measurements are used in the study.

Nyckelord Keyword

(6)

(7)

Abstract

In several projects within ABB there is a need of state and parameter estimation for nonlinear dynamic systems. One example is a project investigating optimisa-tion of gas turbine operaoptimisa-tion. In a gas turbine there are several parameters and states which are not measured, but are crucial for the performance. Such pa-rameters are polytropic efficiencies in compressor and turbine stages, cooling mass flows, friction coefficients and temperatures. Different methods are being tested to solve this problem of system identification or parameter estimation. This thesis describes the implementation of such a method and compares it with previously implemented identification methods. The comparison is carried out in the context of parameter estimation in gas turbine models, a dynamic load model used in power systems as well as models of other dynamic systems. Both simulated and real plant measurements are used in the study.

Keywords: Extended Kalman filter, Parameter estimation, Recursive prediction error method

(8)

(9)

Acknowledgments

First of all I would like to express my gratitude to Geir Hovland and Mats Larsson at ABB Corporate Research in D¨attwil, Switzerland, both for the work related as well as the mountain related guidance.

My examiner, Svante Gunnarsson kindly put me in contact with ABB in the first place.

All the members of the Automation and Control and Utility Solutions groups as well as the other practicants and thesis students at ABB who made my stay in Switzerland a great one.

My family has always supported me on my adventures, also during this Swiss one.

(10)

(11)

Notation

Symbols

h Alarm level in CUSUM test

K(t) Kalman gain

m Number of measurements

n Number of states

P (t) Parameter covariance matrix

Q(t) Process noise covariance matrix, EKF tuning matrix

R(t) Measurement noise covariance matrix, EKF tuning matrix

s Number of parameters

v(t) Process noise at time t

V (t, θ) Cost function at time t

w(t) Measurement noise at time t

x(t) State vector at time t

ˆ

x(tk|tk−1) Estimate of state vector at time tk given data up to time tk−1

y(t) Measurement vector at time t

ˆ

y(tk|tk−1; θ) Prediction at time t given data up to time tk−1 and

parameter vector θ

ε(t, θ) Prediction error vector at time t

λ Forgetting factor

ν Test level in CUSUM test

θ(t) Parameter vector at time t

ψ(t, θ) Gradient of prediction with respect to parameter vector θ

Operators and functions

E Expectation value

∂f

∂x Partial derivative of matrix f with respect to matrix x

arg minθV (tk, θ) Minimising argument θ to function V (t, θ) dx

dθ Sensitivity matrix

⊗ Kronecker product

svd Singular value decomposition

(12)

Abbreviations

AEKF Augmented Extended Kalman Filter

CUSUM Cumulative Sum

EKF Extended Kalman Filter

FACTS Flexible AC Transmission Systems

ODE Ordinary Differential Equation

PEPP Parameter Estimation for Power Plants

RPEM Recursive Prediction Error Method

(13)

viii Contents 6 Experiments 21 6.1 Damped oscillator . . . 21 6.1.1 Model . . . 21 6.1.2 Estimation . . . 21 6.1.3 Discussion . . . 22 6.2 Compressor . . . 25 6.2.1 Model . . . 25 6.2.2 Estimation . . . 26 6.2.3 Discussion . . . 27 6.3 Heat exchanger . . . 30 6.3.1 Model . . . 30 6.3.2 Estimation . . . 32 6.3.3 Discussion . . . 32 6.4 Power system . . . 35 6.4.1 Model . . . 35 6.4.2 Estimation . . . 37 6.4.3 Discussion . . . 48 7 Conclusions 49 7.1 Summary of the results . . . 49

A Algorithms 53 A.1 State-Augmented Extended Kalman Filter . . . 53

(15)

Chapter 1

Introduction

1.1 Background

For the purpose of monitoring and diagnostics of gas turbines there are several parameters of interest that cannot be measured directly in the plants. At ABB parameterised turbine models have been derived which are based on fundamental physical equations. In an ABB project, Parameter Estimation in Power Plants (PEPP), it is investigated how these models can be used together with available measurements from the plants to provide tools for decision making in the plant operation. For instance, the compressor in the gas turbine gets dirty due to particles in the air being compressed. This has the effect that the efficiency decreases. The compressor then has to be washed, but this means that the turbine has to be taken out of operation for a period of time, with a production loss as consequence. The decision when it’s most economic to stop and wash could be made more rational with the monitoring of certain parameters in the compressor. The concept of parameter monitoring also applies to other parts of the turbine plant, such as the heat exchanger used to extract the heat energy in the turbine’s exhaust gases.

In the PEPP project, one of the objectives was, besides applying the estimation methods on the gas turbine model, to develop a toolbox easy to use for state and parameter estimation on arbitrary dynamic systems, since ABB operates in many fields within power and automation. Several toolboxes for parameter estimation exist, for example in the Matlab and gPROMS environments. To our knowledge, however, there are no toolboxes available for tracking time-varying non-linear pa-rameters. The parameters of interest are usually parameters in a physical model (eg. heat and mass balances). Since the toolbox is intended for monitoring pur-poses, the key feature is the ability to track both slow (degradation) and fast (failures) changes in physical meaningful parameters. A black-box parameter es-timation toolbox would not satisfy these requirements. In addition, the existing parameter estimation toolboxes are often cumbersome to integrate into ABB prod-ucts. For these reasons, it was decided to develop a new toolbox from the bottom up.

(16)

2 Introduction When the interest to investigate the use of the toolbox in a power network application arose, it was decided to include this study as a part of the thesis, since it not only would be in line with the initial objective, but could also become a contribution. A brief background of this application follows:

In order to optimise transmission in power networks with respect to cost and safety, it is of interest to monitor the load characteristics at certain points in the network. The relation between load and voltage defines limits to the transferable power at a given working point. When a disturbance occurs in the system, for instance a disconnection of a power line, the new set point should be predicted, so that proper action can be taken in time to keep the system in a safe state. If this is not done correctly, the system could end up either at a stable working point with poor utilisation of the network or, in the other extreme, at an unstable point with a (voltage) collapse as a possible consequence.

1.2 Problem formulation

Given is a model of a dynamic system and data collected from the plant. The model contains parameters to be determined by use of the measured data. The parameters can be constant or time varying. This problem is known as system identification, and more specifically grey-box identification since structural knowledge about the physical system is used when building the model. There exist different approaches to solve this problem, and this thesis will compare two of these: The Extended Kalman Filter (EKF) and the Recursive Prediction Error Method (RPEM).

1.3 Objective

The goal is to implement the RPEM and compare it with the EKF. The comparison will be made regarding the estimation of parameters in test models as well as models used in current ABB projects; Parameter estimation for power plants and Wide area control of FACTS(Flexible AC Transmission Systems).

1.4 Limitations

Apart from the simple example system in chapter 2, the system models used to evaluate the different methods for the parameter estimation will be considered given and no discussion of the background or underlying physics will be made.

1.5 Thesis outline

First, in chapter 2 the concept of grey-box modelling is introduced, followed by a theory section (chapters 3 and 4) where the different filters are described. Im-plementation issues are discussed in chapter 5 and results from simulations are examined in chapter 6.

(17)

Chapter 2

The parameterised model

For the analysis and control of a system it is often of interest to have a mathematical model that describes its behaviour. In any dictionary one can find definitions of systems in the technical sense: A regularly interacting or interdependent group of items forming a unified whole.1 _{If the interaction of these items, and thereby the} system, is governed by laws of physics, for example thermodynamics or mechanics, this knowledge may provide us with mathematical equations to form the model. However, if the system is of a complex nature it may be difficult or impossible to make a detailed model, and one will have to regard the system as a black box where the choice of model structure is based on other approaches. If the insight into the system is complete and the model is a perfect replica of the actual system, it is called a white-box model, and in between the black and the white we find the idea of grey-box modelling. As the shade suggests, this is a model where some, yet not all, of the systems’s structure is known. The unknown part of the model may be expressed in terms of a number of unknown parameters, and different values of these parameters may be tested to complete the model, to make the grey box as white as possible. The goal is to make the model reproduce the behaviour of the true system. This is known as system identification or parameter estimation, to adjust the parameters in the model so that the model reproduces the behaviour of the true system. The same principle applies to the case where the parameters are not constant, in which case the task of the identification procedure is to track the changes of the parameters.

The concept of grey-box modelling will here be illustrated with a simple exam-ple. Consider the damped oscillator in figure 2.1.

Newton’s second law gives us

m¨x = −kx − c ˙x + f (2.1) where x is defined in a coordinate system where x = 0 at the spring equilibrium.

1_{Encyclopaedia Britannica}

(18)

4 The parameterised model

f

m

c

k

Figure 2.1. Damped oscillator.

We define natural frequency and damping ratio as

ω =pk/m and ζ = c

2√km (2.2)

respectively, and write 2.1 on state space form as ˙x = · 0 1 −ω2 _−2ζω ¸ x + · 0 ω2 ¸ f (2.3) y = x (2.4)

The problem we are facing is to determine the values of ω and ζ. At our disposal we have measurements of u and y (possibly corrupted by noise), and we have the mathematical model above, describing the system dynamics. The parameters that are discussed here will be referred to as system parameters.

We will here consider two approaches to solve this problem, by Kalman filtering and by the use of a recursive prediction error method. These two methods have a lot in common, noted in [13] where discrete-time systems were studied, in [14] as well as in [2].

As in [2], from where the algorithms used for implementation have been taken, we will here consider continuous-time dynamic systems with discrete-time mea-surements. A general representation of the state space equations can be written:

˙x(t) = f (x(t), θ(t), u(t)) + v(t) (2.5)

y(tk) = h(x(tk), θ(tk), u(tk)) + w(tk) (2.6)

Here x(t) is the state vector, y(t) is the vector of measurements, u(t) is the system input and v(t), w(tk) are the process and measurement noises respectively, with

(19)

Chapter 3

Kalman filtering

In 1809, as a result of studying the orbits of celestial bodies, Gauss ([4]) wrote: ”. . . since all our measurements and observations are nothing more than approxi-mations to the truth, the same must be true of all calculations resting upon them, and the highest aim of all computations made concerning concrete phenomena must be to approximate, as nearly as practicable, to the truth. But this can be accom-plished in no other way than by a suitable combination of more observations than the number absolutely requisite for the determination of the unknown quantities. This problem can only be properly undertaken when an approximate knowledge of the orbit has been already attained, which is afterwards to be corrected so as to satisfy all the observations in the most accurate manner possible.”

These ideas also constitute the rationale behind the Kalman filter, which was derived on discrete-time form for linear dynamical systems by Kalman in 1960 ([9]). The year after, the continuous-time version was presented by Kalman and Bucy ([10]). We will here give a description of the mixed version of the discrete and continuous filters that is required when dealing with continuous dynamics and discrete measurements.

As a starting point, we will consider the system (2.5-2.6), neglecting for the moment the parameter dependency and assuming that the functions f and h are linear. The general system can then be written

˙x(t) = Ax(t) + Bu(t) + v(t) (3.1)

y(t) = Cx(t) + w(t) (3.2) where A, B and C are matrices of suitable dimensions. The process and measure-ment noise vectors v and w have zero means and variances given by

Q = Ev(t)vT_{(t) and R = Ew(t)w}T_(t) _(3.3)

Given that these noises have gaussian distributions, the best (gives the lowest variance of the prediction error) estimator of the states x(t) given the observations

y(t) is the Kalman filter, which we will describe here.

(20)

6 Kalman filtering The prediction of the state x(tk+1|tk) given the observations y(tk) is obtained

by integrating the process function

˙ˆx(t) = Aˆx(t) + Bu(tk) (3.4)

from tk to tk+1. We denote the covariance of the prediction error

P (t|t) = E£¡x(t) − ˆx(t|t)¢¡(x(t) − ˆx(t|t)¢T¤. (3.5)

P (tk+1|tk) is obtained by integrating

˙

P = AP + P AT _{+ Q} _(3.6)

over the same time interval as above. This is often referred to as the time update or prediction step in the algorithm. As a new measurement is taken into account, the prediction and its corresponding covariance are updated in the measurement update or correction step:

ˆ

x(tk+1|tk+1) = ˆx(tk+1|tk) + K(tk+1)(y(tk+1) − C ˆx(tk+1|tk)) (3.7)

where K is the Kalman gain, given by

K(tk+1) = P (tk+1|tk)CT[CP (tk+1|tk)CT) + R)]−1 (3.8)

and the covariance is updated as follows:

P (tk+1|tk+1) = P (tk+1|tk) − P (tk+1|tk)CT ×

×[CP (tk+1|tk)CT + R]−1CP (tk+1|tk) (3.9)

This algorithm needs to be initiated with values for the state vector x and the covariance matrix P . The choice of these is discussed in 5.6.4. Then the quantities ˆ

x and P can be updated as new measurements become available.

3.1 The extended Kalman filter

When the functions f and h are nonlinear, an often used approach is to linearise the functions, thus achieving a linear system allowing the ordinary Kalman filter to be used. If we make a first order Taylor expansion of f and h in the neighborhood of the latest available state estimate ˆx and take

A = ∂f ∂x ¯ ¯ ¯ x=ˆx and C = ∂h ∂x ¯ ¯ ¯ x=ˆx (3.10)

the recursive algorithm (3.4-3.9) can be used, updating the linearisations (3.10) at every time step using the latest available state estimates. This is the well-known extended Kalman filter (EKF). The implemented version doesn’t use the linearisations in the state propagation (3.4), but uses the non-linear function (2.5) directly.

(21)

3.2 State augmentation 7

3.2 State augmentation

The concept of state augmentation is to include in the state vector unknown system parameters to be estimated. We thereby say there is no formal difference between a state and a parameter. Let the parameters be stacked in the vector θ. We write

¯ x(t) = · x(t) θ(t) ¸ (3.11) The system parameters are modelled as stochastic variables, varying as

˙θ(t) = vθ(t) (3.12)

where vθ is zero mean gaussian noise with covariance Qθ = Evθ(t)vθ(τ )Tδ(t − τ ).

This gives us the possibility to use knowledge about the nature of the parameters, if they are constant or varying, which can be used when tuning the corresponding entries in the augmented process noise matrix ¯Q:

¯ Q = · Qx 0 0 Qθ ¸ (3.13) In this way, augmenting the state vector with the parameters, including (3.12) in the set of process equations and extending the process noise covariance matrix, the Kalman filter (or the extended Kalman filter) can be used to simultaneously track the states and the (time varying) parameters. Note that a system linear in the states often becomes nonlinear using this approach. The convergence properties of this method applied on discrete-time systems was given in [13]. The implemented algorithm of the state augmented extended Kalman filter has been taken from [2] and is given in A.1 for reference.

(22)

(23)

Chapter 4

The recursive prediction

error method

The theory concerning system identification in general is thoroughly described in [14]. The specific part we are interested in is referred to as the Recursive Prediction Error Method (RPEM) and is more specifically covered in [15]. Applied to state space models, this method is closely related to the EKF as pointed out in [13] and [2]. Here we will give a brief description of the RPEM in general and in section 4.1 point out the major differences to the EKF.

This method is as the name suggests based on the prediction error, which is formed by subtracting the predicted output from the measured. Given the past data, i.e., inputs and measurements, the prediction at a certain point in time is obtained from a parameterised model and the algorithm adjusts the parameters in order to minimise the prediction error, thus tuning the model to reproduce the true system’s behaviour.

We take

ε(tk, θ) = y(tk) − ˆy(tk|tk−1; θ) (4.1)

as the difference between the prediction ˆy(tk|tk−1; θ) and the system output y(tk).

We sum these errors over time to measure how well the model describes the true system in a weighted least squares criterion

V (tk, θ) = 1 2 k X i=1 λk−i_εT_(t i, θ)Λ(ti)−1ε(ti, θ) (4.2)

where λ is a forgetting factor used to assign smaller weight to older measurements, further discussed in section 5.6.2. Λ(ti) is a sequence of weighting matrices. At

time tk we choose the parameter vector that minimises this criterion, that is

ˆ

θtk= arg min

θ

V (tk, θ) (4.3)

(24)

10 The recursive prediction error method Using a Gauss-Newton search scheme to find the minimum of this function, includ-ing one more data point in every search iteration gives the followinclud-ing algorithm:

ε(tk) = y(tk) − ˆy(tk) (4.4)

ˆ

θ(tk) = θ(tˆ k−1) + γ(tk)R−1(tk)ψ(tk)Λ(tk)−1ε(tk) (4.5)

R(tk) = R(tk−1) + γ(tk)[ψ(tk)Λ(tk)−1ψT(tk) − R(tk−1)] (4.6)

Here γ(tk) is the adaptation gain. A common choice of forgetting profile is to let

λ(tk) be constant (= λ) as in 4.2. λ and γ are related as follows:

γ(tk) =

γ(tk−1)

γ(tk−1) + λ

(4.7) Equations (4.4) - (4.4) are usually implemented in a different form to avoid the inversion of the matrix R in 4.5. Using the matrix inversion lemma the algorithm takes the following form [14]:

ε(tk) = y(tk) − ˆy(tk) (4.8) ˆ θ(tk) = ˆθ(tk−1) + L(tk)ε(tk) (4.9) L(tk) = P (tk−1)ψ(tk)[λ(tk)Λ(tk) + ψT(tk)P (tk−1)ψ(tk)]−1 (4.10) P (tk) = P (tk−1) − P (tk−1)ψ(tk)[λ(tk)Λ(tk) + ψT(tk)P (tk−1)ψ(tk)]−1ψT(tk)P (tk−1) λ(tk) (4.11) Apart from the parameters, the algorithm also carries their covariance matrix P . Moreover, the gradient of the prediction with respect to the parameter vector,

ψ(tk, θ) = dˆy(tk, θ)

dθ (4.12)

is needed. All we need to provide to be able to use this algorithm is thus a param-eterised predictor and its gradient with respect to the parameters. In the following we will discuss this further.

4.1 Gradients

The basic idea is now to use the Extended Kalman Filter discussed in section 3.1 to track the (parameterised) system’s states. We also parameterise the EKF tuning matrices Q and R. Instead of augmenting the parameters as in section 3.2, we use the algorithm (4.8) - (4.11) with the EKF as the predictor, thus parameterised both through the parameters in the system itself and the parameters in the tuning matrices. Since the prediction is obtained from (2.6), its gradient has the following form dˆy(tk, θ) dˆθ = ∂h ∂ ˆx dˆx(tk) dˆθ + ∂h ∂ ˆθ (4.13)

(25)

4.2 The EKF vs. the RPEM 11 The gradient of the states with respect to the parameters, i.e. dˆx

dˆθ is referred to

as the sensitivity in [2] and its propagation and update equations can be obtained by differentiating the corresponding equations for the states themselves. We will here only account for the derivation of the filter equations up to the point where the interesting connection to the EKF can be seen. For a full derivation of the algorithm, please refer to chapter 4.4 in [2].

We take the state propagation equation (2.5) that is used to propagate the state estimates and differentiate with respect to θ. Using the chain rule this gives the sensitivity propagation d ˙ˆx(t) dˆθ = ∂f ∂x dˆx(t) dˆθ + ∂f ∂θ. (4.14)

For the sensitivity update we begin with the state update equation (3.7). With

ε(tk+1) = y(tk+1) − C ˆx(tk) we differentiate with respect to θ and get

dˆx(tk+1) dˆθ = dˆx(tk) dˆθ − d dθ(K(tk+1)ε(tk+1)) (4.15)

and by using the product rule and noting that the gradient of the prediction error

ε is the negative gradient of the prediction ˆy, we get dˆx(tk+1) dθ = dˆx(tk) dˆθ − K(tk+1) dˆy(tk+1) dˆθ + dK(tk+1) dˆθ ε(tk+1). (4.16)

It is the last term in the update equation, dK(tk+1)

dˆθ , that causes the complexity in

this version of the RPEM. Omitting this term in the continued derivation of the filter gives an algorithm that strongly resembles the EKF ([13], [2]). The differences between the two algorithms will be further discussed in section 4.2 below. The implemented RPEM is taken from [2] and is given in A.2 for reference.

4.2 The EKF vs. the RPEM

The close relationship between these two methods was investigated for discrete-time linear systems by Ljung ([13]) who showed that the missing term dK

dˆθ is what

leads to the convergence problems of the EKF. Ljung started with the EKF and suggested modifications to overcome this problem: To include the term and accept the increased algorithm complexity or to directly parameterise the gain matrix K to make the computations of dK

dˆθ less demanding, yet keeping the good asymptotic

convergence properties of the stochastic descent method.

Also in the thesis of Bohn ([2]) this relationship has been examined, here in the light of the matrix differential calculus applied to write the algorithms on vectorised form. The starting point was the general RPEM and then both the version of the filter with direct parametrisation of the Kalman gain matrix discussed in [13] as well as the full RPEM was derived using his notation.

We can thus see the EKF as the RPEM applied to a state space model with a simplification in the gradient calculation. The simplification neglects the cou-pling between the parameters and the Kalman gain which gives rise to convergence

(26)

12 The recursive prediction error method problems and also makes it impossible to parameterise the EKF tuning matrices as discussed below.

The performance of the EKF strongly depends on the tuning of the covariance matrices Q and R. In the case of a linear system with gaussian noise the optimal filter is achieved when these matrices correspond to the covariance matrices of the noise affecting the process and the measurements. Otherwise these matrices are to be seen as tuning parameters that are adjusted to make the algorithm work in a satisfactory way. In this version of the RPEM the elements in the mentioned matrices are seen as filter parameters to be adapted by the filter itself (along with the system parameters in (2.5-2.6).

The calculation effort needed to compute the term dK

dˆθ in (4.16) i.e. the

deriva-tive of the Kalman gain matrix with respect to the parameters depends on the way the parameters affect the gain. In the Kalman filter this dependence is rather complex, and practically all the filter prediction and update equations (3.4-3.9) have to be differentiated with respect to θ to obtain dK

dˆθ.

To be able to incorporate prior knowledge about the nature of the (time vary-ing) system parameters it is of interest to have the possibility to parameterise the matrices Q and R. We have therefore chosen to implement the version where the filter parameters enter the Kalman filter through the mentioned tuning matrices, instead of directly parameterising the gain. This version of the RPEM is the fil-ter that is compared with the state augmented EKF on the considered systems in chapter 6. The implementation of the RPEM is taken from Bohn ([2]), who calls it an Adaptive Extended Kalman Filter (AEKF). Both the RPEM and the EKF algorithms are given in appendix A.

(27)

Chapter 5

Implementation

Here we will discuss practical issues that has been considered during the imple-mentation of the algorithms in question. We will also give an overview how to tune the filters to alter their behaviour.

5.1 The central processing unit and its efforts

The fact that the recursive prediction error method is computationally demanding was suggested in [13]:. . . will require an amount of computing that may be forbid-ding for higher order systems.”. Since the time of this statement the computer power has increased by a factor of approximately 2(2003−1979)∗12/18 _{= 65536, so} one would think that this no longer presents a problem.1

5.1.1 The EKF

The main computational burden in the EKF algorithm is the solving of the Riccati equation. This equation grows like n2_{, where n is the number of states, and it has} to be numerically integrated between time steps. Measures to be taken can be to remove all the redundant entries in the covariance matrix before the propagation and make sure that the method used for the integration is efficient.

5.1.2 The RPEM

The RPEM suffers from even heavier computations than the EKF. This is due to that first, the second derivatives, or the Hessian, of the system are needed and secondly, the propagation of the sensitivity needs to be numerically integrated. The second derivatives are not a big problem as long as function calls are cheap. The

1_{Gordon Moore stated in 1964 that computational power at a given cost would double every}

12 months, which later has been adjusted to 18 months. On the other hand, the somewhat ironic counterpart to Moore’s law, Gate’s law, states that computer software gets half as fast during that same time period.

(28)

14 Implementation numerical evaluation of the derivatives could of course be replaced by analytical derivatives, but when the functions get complicated, the numerical approach is more appealing in a setting where the user should be required to provide the model only, not its derivatives. The time consuming part is the sensitivity propagation equation, which grows like n + (1 + n)(n + s) where s is the number of parameters. The need for efficient ODE solvers is thus even more important for this filter. Here we can also take away redundant entries in the matrices, but this will only reduce the number by approximately half and not alleviate the ”curse” of dimension.

5.2 Programming environment

The algorithms have been implemented in Matlab and in C++. The transfer to the latter was made when the estimation of the larger systems became too time-consuming for the use of Matlab to be efficient. For this purpose the freely downloadable matrix library Newmat ([3]), has been used which incorporates most of the basic matrix and vector algebraic functions. Unfortunately, the library does not provide the use of sparse matrices, which would have been preferable when implementing the algorithms considered. For the propagation (numerical integration) of states and covariance/sensitivity matrices we have used a fourth order Runge-Kutta solver. We have also tested other, more complex solvers but for the simulation examples in this thesis, no performance gain has been seen compared to the mentioned Runge-Kutta solver. The choice between the different solvers is to be made by the user regarding the nature of the differential equations in the system to be estimated. The solvers listed in table 5.1 are available.

Euler1 Forward Euler integration

Rungekutta4 4th order Runge-Kutta integration

AdaptiveRungeKutta Adapative Runge-Kutta from the Numerical Recipies

in C ([1])

AdaptiveRosenbrock Adaptive Rosenbrock (Stiff solver), also from [1]

LSODE Livermore Solver for Ordinary Differential Equations

(Non-stiff/stiff solver) ([7])

Table 5.1. Tested solvers for numerical integration. The solvers are listed according to their complexity with Euler1 being the simplest.

5.3 Matrix differential calculus

The main contribution from the thesis of Bohn[2] is the application of matrix differential calculus that makes it possible to write all the equations on a completely vectorised form. Matlab is known to be slow when it comes to for-loops and

(29)

5.4 Numerical derivatives 15 therefore Bohn’s implementation with Kronecker products combined with sparse matrices is appealing. Here some of the basics are mentioned to facilitate the reading of the notation in the algorithms. For a complete reference please refer to the mentioned thesis by Bohn.

5.3.1 The Kronecker product

Let A and B be matrices of dimension m by n and p by q respectively. The Kronecker product A ⊗ B is the mp by nq matrix:

A ⊗ B =    A11B . . . A1nB .. . . .. ... Am1B . . . AmnB    (5.1)

5.3.2 Derivatives

Let f (x) be a m by n matrix function and x its p by q matrix argument. Then the derivative df (x)_dx is defined as the np by mq matrix:

df (x) dx = d dx⊗ f (x) =    d dx11 . . . d dx1n .. . . .. ... d dxm1 . . . d dxmn    ⊗ f (x) = (5.2) =     df (x) dx11 . . . df (x) dx1n .. . . .. ... df (x) dxm1 . . . df (x) dxmn    

For the special case where f and x are vectors, _dxdfT becomes the Jacobian matrix

and d_dx2f (x)T_dx is the Hessian matrix.

5.4 Numerical derivatives

The definitions of the first and second derivatives have been taken from Numerical Recipes ([1]), and are defined as the limits when h → 0 of

∂f ∂x =

f (x + h) − f (x − h)

2h (5.3)

for the first derivatives, and

∂f ∂x∂y =

1

4h2[[f (x + h, y + h) − f (x + h, y − h)] −

(30)

16 Implementation for the second. The choice of h in the actual computation has been made according to the discussion in [1] considering machine accuracy.

As different configurations of the derivatives are used in the RPEM, the deriva-tive _dxdxdf is first computed, then restacked to form the needed _dxdfTT_dx and _dxTdf_dxT.

The same goes for the mixed derivatives of _dxdpdf . The first derivatives are simpler to handle, as (_dxdfT)T =

dfT

dx.

5.5 Identifiability

The quality of the estimates is defined by the certainty with which they are deter-mined. How do we know that the estimated parameters are unique? Even if they converge to certain values, it is of interest to know if these values are the only solu-tion or just one among many. A way of monitoring this property of the estimates is to evaluate the rank of the jacobian of the vector function containing 2.5 and 2.6. We stack these and differentiate with respect to the states and parameters and use the following test:

min svd(d[f (t, θ); h(t, θ)]

d[x, θ] )|(x,θ)=(ˆx,ˆθ)< ² (5.5)

and check this for some suitable value of ² > 0. Frequent failures of this test indicate that θ and x can not be assessed from the available measurements. This method is further described in [8]. In this implementation the indicators are computed at even time intervals. A typical execution of the estimation algorithm gives a printout as the one displayed in 5.2, where ”Trust1” is the criterion discussed above while ”Trust2” is a version where the derivation in 5.5 only is made with respect to the parameters.

5.6 Tuning of the filters

The filter performance depends on the choices of tuning parameters such as initial values, weighting matrices etc., different choices to be made by the user. The time spent on tuning the filters must of course stand in proportion to the expected performance gain. Moreover, if a filter is extremely difficult to tune to a satisfactory behaviour, there is also the question how robustly it can be implemented in a real-time setting.

If the parameters we are trying to estimate are not constant, we need to have this in mind when tuning the filters. There is always a fundamental tradeoff between a filters tracking capabilities and its robustness to noise and disturbances. In other words, if we want a filter to respond quickly to changes, the estimates will be ”noisy”.

First we will discuss filter specific tuning issues, followed by an overview of filter independent choices. We will also indicate how to adjust the filters to comply with changing parameters.

(31)

5.6 Tuning of the filters 17

5.6.1 The EKF

The main design variables of the EKF are the elements of the matrices Q and

R, which correspond to the covariances of the process and measurement noise

respectively. For a linear system with gaussian noise the true covariance matrices give the best performance, but this is not true for the non-linear case, where these matrices are to be seen as tuning parameters to obtain a satisfactory behaviour of the filter. In addition, wrongly tuned covariance matrices can give biased parameter estimates ([13], [2]).

When looking at equation 3.6, it is clear that large values of Q gives a large

P and in turn a large Kalman gain, K. This means a large step size with which

Estimating Parameters with Filter: AEKF Individual element adaptation of Q diagonal Individual element adaptation of R diagonal Number of filter parameters: 4

Number of estimated system parameters: 1 Size of sensitivity propagation matrix: 2x17

5 % completed. t=2.5 Trust1=868.816 Trust2=0.0294549

10 % completed. t=5 Trust1=877.919 Trust2=0.0293159

100% completed. t=50 Trust1=870.232 Trust2=0.0295367

Estimation Complete

(32)

18 Implementation the parameters2 _{are updated in 3.7, and the filter adjusts the parameters faster} to track changes. Sometimes there is prior knowledge about the nature of the parameters, for example that different parameters change more rapidly than others. This information could and should be used when tuning the Q matrix. The R matrix represents the confidence in the measurements, a large value means that we assume that the measurements are corrupted by noise with a large covariance, and will make the gain K smaller in 3.6.

5.6.2 The RPEM

Theoretically, every element of the Q and R matrices can be parameterised as long as the number of unknowns does not exceed nm, where n is the number of states and m the number of measurements ([2]). However, since the complexity of, and thereby the computation time needed to run the algorithm, heavily depends on the number of parameters, it is of interest to keep this number as small as possible. Therefore, in this implementation the user can choose between the following con-figurations: Every element of the diagonal of a matrix is parameterised, the matrix is parameterised using one parameter times an identity matrix, or the matrix is not parameterised at all. Every combination of the above for each of Q and R is also valid. If there is prior knowledge of some suitable scaling of the diagonals, this could also be used instead of the identity matrix in the one-parameter approach option. The approach used in the experiments has been to first test with fully parameterised diagonals, then lower the number of parameters until a decrease of performance is noticed.

Using a forgetting profile that gives more weight to recent measurements keeps the algorithm focused on the present behaviour while ”forgetting” older data. This enters the RPEM in the tuning parameter λ. This parameter affects the whole algorithm and it is not possible to include knowledge about differences in variation speed between separate parameters.

Bohn [2] used prior knowledge about which parameter that would change and when, using this to apply a lower bound on the corresponding element in the parameter covariance matrix P just before the changes occurred. This method has been implemented with the lower bound as a tuning parameter.

Ljung [14] suggests a method that resembles the Kalman filter state augmen-tation approach, to introduce a model of the parameter changes. The assumption that the parameters vary according to

˙θ = vθ with EvθvTθ = Qθ (5.6)

will have the impact on the RPEM algorithm that Qθ is added in the parameter

covariance update equation (4.11), corresponding to the update equation of P in the implemented algorithm. This has thus the same effect as the bounding discussed

2_{What said about parameters here also apply to the system states, since there is no formal}

dif-ference between states and parameters in the EKF. We choose, however, to refer to the parameters in this context.

(33)

5.6 Tuning of the filters 19 t T_λ λ end λ start λ 0

Figure 5.1. Plot of lambda over time.

above, to keep the parameter covariance matrix from getting too small, but this approach makes it easier for the user since the tuning becomes identical to that of the EKF. Known differences between different parameters are in this way easy to incorporate in the filtering process, and we therefore propose this change to the implemented algorithm, instead of the lower bounding of P discussed earlier.

In [14] as well as in [2] a forgetting factor, λ, is proposed for two purposes: first to make the filter handle varying parameters and secondly, to improve its transient behavior. These approaches are combined in the following way of computing λ at time tk:

λ(tk) = λend− (λend− λinit)e−

7tk

Tλ _(5.7)

in this way, we can specify the start and end value, λinit and λend, and

approx-imately over how many data λt will increase from λinit to λend, Tλ. Keeping

lambda small in the beginning of the data set increases the initial convergence rate. A slightly different version of this is mentioned in [14] where it is motivated: ”The reason is apparently that early information is somewhat misused and should carry a lower weight. . . ”. The value λendcorresponds to the λ mentioned in section

4. The factor (7) in the exponent has been chosen to give Tλ the desired length

scale. The principal appearance of lambda is displayed in Figure 5.1, and typical values [14] for the parameters are λend∈ [0.98 0.995] for a tracking filter, λend= 1

for estimation of constant parameters and λinit∈ [0.95 λend].

5.6.3 The CUSUM test

If a filter is working well, the residuals of the filter should be white noise. This indicates that all information has been extracted from the available data, thus the filter is doing its very best. A way of testing whether a signal has an average that exceeds a certain level is the CUSUM test described in [6], and the algorithm is given here:

(34)

20 Implementation After a change in the parameters, the residuals will no longer be white, and their average will be larger than the level ν. After an alarm, g(t) is set to zero, and appropriate action is taken to give the filter a chance to re-converge under the new conditions. For the EKF, this action may be to change the tuning covariance matrices Q and R in some fashion and for the RPEM to reset the parameter co-variance matrix P to its initial (or some other) value. The extra design parameters that have to be tuned are h and ν.

5.6.4 Initial values

The filters estimate the states of the system as well as system, and, in the RPEM case, filter parameters. The filters also carry the parameter and state variances, and all these quantities need initial values for the first recursion. If typical values for the parameters are known, these should be used. The start-up covariance matrix will affect the step size in the beginning of the estimation and should reflect the certainty of the initial parameter guess. In the RPEM, initial values for the filter parameters in the tuning matrices Q and R must also be given and often these values are very uncertain. This should be reflected in the corresponding elements of P which in most cases should be kept significantly larger than the elements coupled to the system parameters.

(35)

Chapter 6

Experiments

Here the filters have been applied to a number of dynamical systems. For each sys-tem we first give a description of the model, then we account for how the data sets have been generated. The two filters are then applied to estimate the parameters, followed by a discussion of the results.

6.1 Damped oscillator

6.1.1 Model

The estimation algorithms are here applied to the test system described in section 2.

6.1.2 Estimation

The system has been simulated using the parameter values ω = 10 and ζ = 0.1. The noise added during simulation was gaussian with covariance matrices

Q = · 1 0 0 1 ¸ R = · 1 0 0 10 ¸ (6.1) Up to a certain level of process and measurement noise the EKF gives correct estimates of the parameters. When more noise is included in the simulation, the EKF estimates becomes biased as shown in 6.1. The estimation was also carried out for a longer period of time to verify that what we see is not only slow convergence. In figure 6.2 the simulated (noisy) measurements together with the predictions are plotted. The difference of these two signals, the prediction error, should be white noise when the filter is working, since they represent what the model can not predict. This signal is displayed in figure 6.3, which is an indication of good filter behaviour.

Under the exact same conditions the RPEM parameter estimates converge to the correct values as shown in figure 6.4. The plots of measurements and RPEM

(36)

22 Experiments prediction errors have been left out since they look identical to those of the EKF.For this estimation the diagonals of both matrices Q and R were parameterised giving a total of four filter parameters. These parameters are shown in figure 6.5.The time needed for the estimations were measured and not surprisingly the EKF beats the more complicated RPEM on this point. To estimate 100 seconds of system activity, the EKF finished in about 15 seconds, the RPEM needed 120 seconds.

6.1.3 Discussion

The EKF gives slightly biased estimates under certain conditions involving high noise levels, still the prediction errors are white noise. The interesting point is that the RPEM gives correct estimates under the exact same conditions. Even for this small system the higher complexity of the RPEM results in an execution time of approximately eight times that of the EKF.

10 20 30 40 50 60 70 80 90 100 9 9.5 10 10.5 Time (seconds) Parameter Estimates 10 20 30 40 50 60 70 80 90 100 0.05 0.1 0.15 0.2 0.25

Figure 6.1. Biased EKF parameter estimates. (The correct values are 10 and 0.1 respectively.)

(37)

6.1 Damped oscillator 23 20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 −5 0 5 Time (seconds)

Measurements and filtered measurements

20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 −20 −10 0 10 20

Figure 6.2. Measurements (noisy) and EKF predictions (or filtered measurements) of position (top) and velocity of the spring.

0 5 10 15 20 25 30 35 40 45 50 −4 −2 0 2 4 Time (seconds) Prediction Errors 0 5 10 15 20 25 30 35 40 45 50 −15 −10 −5 0 5 10 15

(38)

24 Experiments 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 Time (seconds) Parameter Estimates 0 5 10 15 20 25 30 35 40 45 50 0 0.05 0.1 0.15 0.2 0.25

Figure 6.4. The parameter estimates obtained with the RPEM

0 5 10 15 20 25 30 35 40 45 50

0 0.5 1

Time (seconds)

Estimated covariance matrix elements(Q and/or R)

0 5 10 15 20 25 30 35 40 45 50 0 0.5 1 0 5 10 15 20 25 30 35 40 45 50 0 0.5 1 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6

(39)

6.2 Compressor 25

6.2 Compressor

The compressor that delivers the pressurised air for the combustion in the gas turbine (figure 6.6) gets dirty due to particles in the input air. This causes an efficiency reduction which has to be alleviated through a wash of the compressor. The wash in turn means that the compressor (thus the entire turbine) has to be shut down for a period of time, an event that obviously costs money due to the production loss. On the other hand, low efficiency also costs money and this is the tradeoff: At what level of efficiency is it economically justifiable to make the wash? The PEPP project aims to develop a tool that addresses this optimisation problem, in which an estimate of the compressor efficiency over time is crucial.

Figure 6.6. Image of a gas turbine with the compressor located at the far right.

6.2.1 Model

A parameterised model of the compressor taken from [5] has the state equations ˙ m(t) = A L ¡ p0(t) − p(t) ¢ (6.2) ˙p(t) = a 2 Vp ³ m(t) − kt p p(t) − pamb(t) ´ , (6.3)

(40)

26 Experiments and the following measurement equations

tac(t) = tis(t) − tamb(t) η + tamb(t) (6.4) tit(t) = ff(t)Cg(t) + cptac(t)m(t) ff(t)cp2+ m(t)cp , (6.5) where tis(t) = tamb(t) µ p0(t) pamb(t) ¶κ−1 κ . (6.6)

For further reference to the model please refer to [5]. The quantities in the equations above are listed in table 6.1.

States: m(t) Compressor mass flow

p(t) Plenum pressure

Inputs: pamb(t) Ambient pressure

p0(t) Pressure downstream compressor

tamb(t) Ambient temperature

Cg(t) Fuel heat value

ff(t) Fuel mass flow

Measurements: tac(t) Temperature after compressor

tit(t) Turbine inlet temperature

Compressor specific A Duct area

constants: L Duct length

a Inlet stagnation sonic velocity

Vp Plenum volume

kt Throttle gain

κ Specific heat ratio

η Isentropic efficiency

Physical constants: cp Specific heat for air

cp2 Specific heat for gas

Table 6.1. List of the quantities in the compressor model 6.2-6.6

6.2.2 Estimation

Real plant measurements have been used for this estimation. As discussed above, the parameter we want to estimate is the compressor efficiency, η. The EKF

(41)

6.2 Compressor 27 estimate is plotted in figure 6.7. The RPEM gives a somewhat smoother estimate shown in figure 6.8. When we look at the predictions and measurements we can see that the RPEM (figure 6.10) reproduces one of the measurements better than the EKF (figure 6.9).

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9

Figure 6.7. Estimated efficiency (η) with the EKF

6.2.3 Discussion

Here we have tested the filters on the compressor model where we have estimated the efficiency to be used in the compressor wash optimisation algorithm. The filters show a slight difference regarding the output predictions, but it is hard to say whether this is a matter of tuning or not.

(42)

28 Experiments 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9

Figure 6.8. Estimated efficiency (η) with the RPEM

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 750 800 850 900 950 1000 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 1350 1400 1450 1500 1550 1600 1650

(43)

6.2 Compressor 29 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 780 790 800 810 820 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 1320 1340 1360 1380 1400 1420 1440

(44)

30 Experiments

6.3 Heat exchanger

The exhaust from the gas turbine is led through a heat exchanger to make use of the surplus heat energy. This is done to increase the overall efficiency of the plant. It is here investigated how the filters can be used to detect faults in the heat exchanger by monitoring of certain parameters. A fault is simulated by a change in the heat transfer coefficient between steam and metal, Ksand this example will

be used to study different ways of dealing with parameter changes. Only simulated data are used.

6.3.1 Model

The heat exchanger model is built up from sub models, each sub model representing a pipe. The configuration of a model with two pipes is illustrated in figure 6.11.

Gas Exhaust gas -Q1 Metal Tm1 -Q1 Steam Steam ? Ps1 ? Ts1 6fs1 6 Pg1 Tg1 6 ? fg1 Gas -Q2 Metal Tm2 -Q2 Steam ? ? 6 ? Ps2 ? Ts2 6fs2 6 6 ? 6 Pg2 Tg2 6 ? fg2 Pipe 1 Pipe 2

Figure 6.11. Illustration of heat exchanger model

The dynamics of each pipe are governed by the following equations:

˙ρ(t) = 1 Vs ¡ fi(t) − fo(t) ¢ (6.7) ˙ Ts(t) = 1 ρVs ¡ Tsi(t)fi(t) − Ts(t)fo(t) + Q(t)/cs− ˙ρ(t)VsTs(t) ¢ (6.8) ˙ Tg(t) = 1 ρgVg ³¡ Tgi(t)− Tg(t) ¢ fg(t) − Q/cg ´ . (6.9)

(45)

6.3 Heat exchanger 31 Two states are measured, Tsand Tg, and we also have two measurement equations:

fs(t) = fr s Pi(t) Po(t)− 1 (6.10) Tm(t) = Kgfg(t) 0.6_T g(t) + Ksfs(t)0.8Ts(t) Kgfg(t)0.6+ Ksfs(t)0.8 (6.11) and static relations

Q(t) = Kg ¡ Tgi(t) − Tm(t) ¢ fg(t)0.6 (6.12) Tgo(t) = Tgi(t) − Q(t) fg(t)cg. (6.13)

The variables and constants are listed in table 6.2

States: ρ(t) Steam density

Ts(t) Steam temperature

Tg(t) Gas temperature

Inputs: Tgi(t) Input gas temperature

Tsi(t) Input steam temperature

fg(t) Gas flow

Pi(t) Input steam pressure

fo(t) Output steam flow

Measurements: Tm(t) Metal temperature

fi(t) Input steam flow

Signals in static relations: Tgo(t) Gas temperature after heat exchanger

Q(t) Heat from gas path

Constants: ρg Gas density

Vs Steam volume

Vg Gas volume

fr Friction coefficient

cs Specific heat for steam

cg Specific heat for gas

Po Output steam pressure

Kg Heat transfer coefficient gas to metal

Ks Heat transfer coefficient metal to steam

(46)

32 Experiments

6.3.2 Estimation

Only simulated data have been used in this estimation. Measurement noise has been added. Here we estimate the heat transfer coefficients Ks and Kg in two

pipes. At time 1000 a change in Ksin the first pipe is simulated, corresponding to

a sudden fault in the pipe. The parameter changes from a value of 8000 to 6500. Two different tunings of the EKF have been tested to point out the tradeoff between tracking capability and noise sensitivity. The parameter estimates from a slow EKF are shown in figure 6.12. The filter tracks the change correctly. The result from a faster filter is displayed in figure 6.13. Here the filter reacts more rapidly to the change, at the cost of a higher noise level in the parameter estimates. The RPEM has been tested using two approaches for the tracking of the changing parameter. In figure 6.14 bounding of the matrix P has been used and this results in a behaviour similar to the faster EKF. The filter tracks the change quite well, but the estimates are noisy. The CUSUM test has been tried and the resulting parameter estimates can be found in figure 6.15. Here the estimates do converge nicely yet the filter manages to track the change in the first parameter. On the other hand we can see that the estimates of the constant parameters are also affected by the changes. In these estimations all parameters have been treated the same way. Knowledge about which parameter that would change has not been used, since it would be unrealistic to assume that one would know in which of the pipes the fault would occur.

The average execution times for these estimations were 30 seconds for the EKF and 890 seconds for the RPEM.

6.3.3 Discussion

The fundamental tradeoff between the tracking capability and the noise robust-ness has been displayed. Both filters give correct estimates of the parameters and handle the case where parameters change. Since the RPEM is designed to esti-mate constant parameters, measures have to be taken when parameters change. Two approaches have been evaluated, to apply a bound to the matrix Ppp and the

CUSUM test. The CUSUM test turned out to be rather difficult to tune, but when it is working satisfactorily it is clearly an improvement to the algorithm. The es-timation with the RPEM takes approximately 30 times longer than with the EKF for this system size.

(47)

6.3 Heat exchanger 33 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6000 7000 8000 0 200 400 600 800 1000 1200 1400 1600 1800 2000 7500 8000 8500 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6500 7000 7500 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6500 7000 7500

Figure 6.12. EKF parameter estimates

0 200 400 600 800 1000 1200 1400 1600 1800 2000 6000 7000 8000 0 200 400 600 800 1000 1200 1400 1600 1800 2000 7500 8000 8500 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6500 7000 7500 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6500 7000 7500

(48)

34 Experiments 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6000 7000 8000 0 200 400 600 800 1000 1200 1400 1600 1800 2000 7500 8000 8500 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6500 7000 7500 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6500 7000 7500

Figure 6.14. RPEM parameter estimates with bounding of P

0 200 400 600 800 1000 1200 1400 1600 1800 2000 6000 7000 8000 Parameter Estimates 0 200 400 600 800 1000 1200 1400 1600 1800 2000 7500 8000 8500 9000 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6500 7000 7500 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6500 7000 7500

(49)

6.4 Power system 35

6.4 Power system

For protection, control and analysis of power systems it is of interest to be able to estimate the characteristics of the connected load at certain points in the net-work. To be able to operate closer to the limits thus transferring as much power as possible is interesting from an economical point of view, especially with the present deregulation of the power supply utilities that presents new demands on cost efficiency. The structure and voltage levels of one part of the power system in southern Sweden is shown in figure 6.16. The points we are interested in are the load buses that are indicated by an arrow with the text ”load”. The active and reactive load demands acting on such a bus are considered to be dynamically related to the bus voltage. Hill and Karlsson [11] have proposed an aggregated load model that describes this load-voltage relation in the bus. Here we will make a feasibility test to see if the parameters in their model can be estimated with the studied filters, and if so, compare the filters in accordance with the objective of the thesis.

Figure 6.16. The model structure of a part of the power network in the south of Sweden. (The abbreviations stand for Tomelilla, J¨arrestad and ¨Ostra Tommarp.)

6.4.1 Model

Since the load consists of a great number of different devices such as lighting, heating and motors and the measurements are affected by tap changers and other control devices it is virtually impossible to construct a model based on physical laws due to the complexity and size of the system. The dynamical model considered in [11] is a combination of an aggregated static power dependency of voltage that describes the response in active or reactive power after a change in the voltage:

P (t) = P0 ³ V (t)

V0 ´α

(6.14) and first order dynamics that describe the observed power recovery. For further reference to the origin of the model, please refer to [11]. The model has the following

(50)

36 Experiments structure: ˙xp(t) = P0 ·³_{V (t)} V0 ´αs −³ V (t) V0 ´αt¸ −xp(t) Tp (6.15) ˙xq(t) = Q0 ·³_{V (t)} V0 ´βs −³ V (t) V0 ´βt¸ −xq(t) Tq (6.16) Pl(t) = xp(t) Tp + P0 ³ V (t) V0 ´αt (6.17) Ql(t) = xq(t) Tq + Q0 ³ V (t) V0 ´βt (6.18) Note that the equations for active and reactive power are identical and decou-pled. Unlike the active, the reactive nominal load Q0is not always positive, and the situation where it is zero, or close to zero can give rise to problems as mentioned in [16]. The quantities in the model are listed in table 6.3.

States: xp(t) Active power recovery

xq(t) Reactive power recovery

Inputs: V (t) Supply voltage

Measurements: Pl(t) Active power consumtion

Ql(t) Reactive power consumption

Parameters: V0 Pre-fault supply voltage

P0 Active power consumption at pre-fault voltage

Q0 Reactive power consumption at pre-fault voltage

αs Steady state active load-voltage dependance

αt Transient active load-voltage dependance

βs Steady state reactive load-voltage dependance

βt Transient reactive load-voltage dependance

Tp Active load recovery time constant

Tq Reactive load recovery time constant

Table 6.3. List of quantities in the load model

This model has been identified in field tests in the power network in southern Sweden in [11] and an approach with a linearised model and automated identifica-tion has been developed in [16]. These identificaidentifica-tion procedures are both based on analysing a step response in the load after a step has been applied in the voltage. Then V0, P0and Q0 are chosen as the pre-step values of V , Pland Qlrespectively.

In the first approach the parameters αt was then computed using the value of Pl

(51)

6.4 Power system 37 using a least-squares method (The same goes for the corresponding reactive part of the model). In the second approach the system was first linearised around the working point U0 and then the linear parameters were estimated using a least-squares method and then the original parameters were extracted from the linear estimates through an iterative scheme. The automation in this approach consists of detecting steps of a certain quality in the input, i.e the voltage, and thereafter fitting the model to the step response in the mentioned manner.

The study here aims to investigate if it is possible to use methods that would be continually active, i.e., not event triggered and estimate the parameter using the ever-present small variations in the measured voltage. The reason for estimating the parameter values is to be able to forecast the steady-state level of the loads at a time directly after the step. This will let us predict a dangerous situation in the network, and then use the time before the load recovers to take action to prevent an unwanted situation. It is clear that a step caused by a tripped transmission line is not a suitable starting point for the estimation since that would be the time to have access to the estimated parameters. For on-line use, it is thus imperative that it be possible to estimate the parameters using only minor excitations in the voltage. In [11], steps of a magnitude of 10% of the initial voltage level were used. The automated scheme in [16] triggered the estimation procedure on voltage steps of 1.6% or higher.

6.4.2 Estimation

No complete sets of real data from actual operation have been available for this estimation. We have studied two ways to generate the data, referred to here as test cases. In the first test case a simplified system (6.17) has been used to generate the data in Modelica, where only step signals have been used as system inputs and to simulate parameter changes. The second test case has been a model of the system in figure 6.16 which is a more realistic structure and here we have used real data to some extent: The nominal load levels have been taken from real measurements, but the voltage and load signals are generated by simulation. Further information about this system modelled in Modelica can be found in ([12]). In this study the two filters have shown a very similar behaviour. In an effort to give the testing some structure we have also divided the procedure in three setups where different combinations of the parameters have been estimated. Moreover, since there would be no point in showing two versions of almost identical plots from the different simulations we have chosen just to display results from one of the filters for every setup.

First setup

According to the law of parsimony1 _{we start by trying the estimation on the first} test case keeping V0, P0and Q0constant set to 1. The other parameters are chosen as typical values listed in table 6.4.

(52)

38 Experiments

Active load αs = 1.28 αt = 1.92 Tp = 173

Reactive load βs = 1.32 βs = 2.92 Tq = 83.5

Table 6.4. Typical values for the load model parameters.

The input used is a sequence of steps (figure 6.19) with a step size of 2%. All the parameters converge to their correct values in this setup (figure 6.18). The variance of the parameters are shown in figure 6.20 where we can see that the values corresponding to the parameters Tpand Tq(the two bottom curves) decrease slower

than the others. This indicates that the data contain less information about these parameters. Note that for the other parameters there is a sudden decrease at the time of the negative flank of the first step in the input. Here we only let the RPEM adjust the elements of the process noise covariance matrix Q while keeping R constant, to keep the number of parameters down. There was no notable difference from estimating both matrices. The estimated diagonal of Q is shown in figure 6.21.

For the second test case the input is shown in figure 6.22 and the parameter estimates are displayed in figure 6.23. Also in this case all the unknown parameters converge nicely. Figure 6.24 shows the estimated Q matrix elements.

All the estimates in this setup are generated with the RPEM, but very similar results (therefore omitted) were obtained using the EKF.

(53)

6.4 Power system 39 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0.5 1 1.5 Time (seconds) P0 Parameter Estimates 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0.5 1 1.5 Q0 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 2 4 αs 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 2 4 αt 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 −5 0 5 βs 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 5 βt 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 149 150 151 Tp 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 95 100 105 Tq

Figure 6.18. Parameter estimates, first setup with constant nominal loads. All param-eters converge to correct values.

0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0.98

1 1.02