Data Driven Fault Tolerant Control: A Subspace Approach

(1)

Data Driven Fault Tolerant

Control: A Subspace Approach

(2)

(3)

CONTROL: A SUBSPACE APPROACH

PROEFSCHRIFT

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus Prof. dr. ir. J.T. Fokkema,

voorzitter van het College voor Promoties,

in het openbaar te verdedigen op

woensdag 11 november 2009 om 10:00 uur

door

Jianfei DONG

Master of Engineering, National University of Singapore

geboren te Shijiazhuang, Hebei, China

(4)

Prof. dr. ir. M. Verhaegen

Samenstelling promotiecommisie:

Rector Magnificus, voorzitter

Prof. dr. ir. M. Verhaegen, Technische Universiteit Delft, promotor Prof. dr. R. Babuˇska, Technische Universiteit Delft

Prof. dr. ir. J. Hellendoorn, Technische Universiteit Delft Prof. dr.-ing. S. X. Ding , Universität Duisburg-Essen Prof. dr. F. Gustafsson, Link öping University Prof. dr. M. Kinnaert, Université Libre de Bruxelles Prof. dr. ir. E. Holweg, SKF & Technische Universiteit Delft Prof. dr. ir. J. A. Mulder, Technische Universiteit Delft, reservelid

This dissertation has been completed in partial fulfillment of the requirements for the graduate study of the Dutch Institute of Systems and Control (DISC).

The work presented in this thesis has been supported by the AI4IA project of the European Marie Curie FP6 Research Training Programme: 514510 and SKF, the Netherlands.

ISBN: 978-90-9024772-4

Copyright c 2009 by Jianfei Dong.

All rights reserved. No part of the material protected by this copyright notice may be re-produced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any information storage and retrieval system, without writ-ten permission from the copyright owner.

(5)

(6)

(7)

Acknowledgments

T

ime flies. As the four-year research in Delft approaching to the end, my long career as a student is most likely going to a full stop. I still remember the inspiration I got during my high-school years from a book telling the stories of some great scientists, like Marie Curie. Years later, when doing my Master’s study in Singapore, I decided to do a PhD in the West, the cradle of the modern science and technology. This dream became true in 2005, when Prof. Michel Verhaegen offered me the opportunity of working on the European Marie Curie project.

My thanks to Michel are certainly not only limited to this opportunity. Dur-ing the four years, Michel devoted a lot of efforts in supervisDur-ing my study. I will not forget the numerous ideas he brought up to me; will not forget the time he spent with me in debugging the simulation program for identifying a wind tur-bine model; and will not forget the encouragement he gave me, when our paper to SafeProcess09 was rejected due to a misunderstanding in a review. Thank you, Michel, for the knowledge that you have conveyed to me, and for the scientific thinking that you have steered me to.

During the first three years, my research was sponsored jointly by the Euro-pean Marie Curie FP6 Research Training Programme under the project number 514510 and SKF Engineering & Research Center, in the Netherlands. In the course of this project, I got the opportunities in cooperating with leading European in-stitutes, and attended a series of training programs on artificial intelligence and its industrial applications. For these, I would like to thank the project coordina-tors, Dr. Slavco Velikov and Dr. Armen Laziev, for their efforts in organizing this project. I would also like to thank Dr. Edward Holweg for arranging my visit to SKF-Airasca, in Italy in 2006. Many thanks also go to Dr. G ¨oran Christiansson and Dr. Adam Reedman for their arrangement and kind help during the assignment on the power-assisted lifting device in SKF. I would also like to thank Dr. Florin Tatar, for his kind offer of riding with him from Delft to SKF in Nieuwegein; and like to thank Dennis de Bot, for his kind help in soldering the cables for the setup in SKF.

Then I would like to thank Dr. Balazs Kulcs´ar, with whom I have been coop-erating closely since December 2007. Special thanks certainly go to Dr. Xiukun Wei, whose help and friendship helped me go through a very difficult period of my time in Delft. Besides, I will not forget the help of Dr. Stoyan Kanev and Dr. Daniele Corona during the first two years of my study in Delft. Without Daniele’s

(8)

friendship, It would have been more difficult for me to get into the community in Europe. I would also like to express my sincere appreciation for Dr. Jan-Willem van Wingerden and Diederick Joosten, for their help in translating the summary into Samenvatting. Thanks will also go to Mw. Kitty Dukker, Mw. Ellen van den Berg-Moor, Dr. Rufus Fraanje, Dr. Xavier Bombois, Dr. Ton van den Boom, Dr. Yucai Zhu (TU Eindhoven), Dr. Redouane Hallouzi, Dr. Rudy Negenborn, Math-ieu Gerard, Stefan Kuiper, Lakshmi Baskar, Edo de Bruijn, and Ivo Houtzager for their valuable discussions and help. I would definitely also like to thank Xiali Wang, Hong Song, Jing Shi, Shu Lin, and all my Chinese friends in Delft, for their help and friendship.

Appreciations must also be expressed to all the members of my PhD defense committee for their time investment in reviewing my thesis and their valuable comments.

Finally, I would like to thank my parents and brother, for their love, care, sup-port, encouragement, wishes, and understanding, through all the 12 years since 1997, when I left home and pursued my studies. They certainly deserve sharing the accomplishment of my PhD study!

Delft, October 2009 Jianfei Dong

(9)

1 C

Introduction

A

controlled process generally consists of actuators, control laws (e.g. software) and sensors. Faults in one of these components normally lead to unacceptable performance of the system. In many cases, the fail-ure itself is not that critical, it is the combination with the control laws or the feedback elements that might lead to an unstable system or unac-ceptable system performance. By detecting the failure and subsequently adapting the control laws, the performance might be adjusted to a level that is still within acceptable limits. This scheme is known as active Fault Tolerant Control (FTC), which usually consists of fault detection and iden-tification (FDI) and controller reconfiguration (CR). Unlike the existing classical model-based approaches, this thesis aims at developing a uni-fied framework of identification, control, fault detection, and persistently exciting input design. These problems are analyzed and developed in a common data equation format, whose dimensionality can be chosen with a limited number of tuning levers. This common modeling format and the limited design parameters empower a systematic treatment and more usability in the whole identification, control design, monitoring, and ex-periment design cycle. The improved usability is a further consequence of the data driven approaches; i.e. to maximize what can be retrieved from the data without the need that user should interfere.

1.1 Introduction to fault tolerant control

Highly automated systems are among the symbols of the modern world. The importance of automation lies in the increased productivity, reliability, and prof-itability, which it brings to an industrial process; and the comfort and ease that it brings to people’s life. Modern technical processes are becoming highly integrated systems of electronic, electrical, mechanical, optical, and perhaps even biological components. These different components exchange information and interact with each other. This, however, also significantly increases the complexity of the system

(14)

at the same time. A failure in any of these components may cause the closed-loop controlled system to behave in an undesirable way.

For an industrial process, this may cause damage to the plant equipment and increase in the wasteful use of raw material and in the downtime of the produc-tion process. According to a survey inNimmo(1995), abnormal situations account annually for at least$10 billion in lost revenue in the U.S. alone. In mechanical

industries, tool failure resulted from wear represents about 20% of machine tool down-time, and negatively impacts the work quality in the context of dimensions, finish, and surface integrity (Liang et al. 2002). In the case of civil applications, fail-ure in a system may jeopardize lives. For instance, the multiple engine separations on the right wing caused the crash of the Boeing 747-200F freighter on 4 October 1992 in Amsterdam, which claimed 43 lives (Maciejowski and Jones 2003).

To be more precise, we shall refer to the following definition of “fault”. Definition 1.1 (Fault (Isermann and Ball´e 1997)) Fault is an unpermitted deviation of at least one characteristic property or parameter of the system from the acceptable/usual/ standard condition.

Since a system usually contains multiple components and has a variety of prop-erties, a fault may occur on any part of the system. As far as a controlled system is concerned, the following figure categorizes the potential faults therein (Blanke et al. 2003).

u Actuators Plant Sensors Actuator faults Plant faults Sensor faults y

Figure 1.1:Types of faults in a controlled system.

To be more precise, the faults in a controlled system can be classified as follows (Blanke et al. 2003).

• Plant faults: Such faults change the dynamical properties of the system. • Sensor faults: The plant properties are not affected, but the sensor readings

have substantial errors.

• Actuator faults: The plant properties are not affected, but the influence of the controller on the plant is interrupted or modified.

From modeling perspective, faults can be classified as follows (Kanev 2004). • Additive faults: Such faults appear as additional terms in the state

equa-tions (see Section1.2of this thesis for the notion of a state equation) of a system. For stochastic models, they result in changes of the mean value of the measured signals only.

(15)

• Multiplicative faults: These faults correspond to changes in the parameters of the state equations or changes in the variance of the stochastic disturbance and noise.

A controlled system needs two functionalities to tolerate these potential faults (Blanke et al. 2003).

• Fault detection and identification (FDI): The existence of a fault has to be detected and the fault identified.

• Control reconfiguration (CR): The controller has to be adapted to the faulty situation so that the overall system continues to satisfy its goal.

With these two units, an FTC system should have the architecture as schematically shown in Fig.1.2(Blanke et al. 2003).

u Controller y faults Diagnosis Controller reconfiguration faults yref Supervision level Execution level Plant

Figure 1.2:The architecture of fault tolerant control.

In what follows, modeling techniques embedding faults into state space mod-els are first discussed with examples. Then the literature on model-based FDI and CR is reviewed. Classical designs based on identified models are introduced in Section1.5. The insight on the role of the data equation in identification, control, and detection motivates the unified design, and its data driven solutions. This chapter ends with the scope, organization, and contributions of this thesis.

1.2 Dynamic models contaminated with faults

1.2.1 LTI state-space models with faults

In modern control theory, LTI systems are usually modeled in the following state-space form,

˙x(t) = Acx(t) + Bcu(t)

(16)

Here x(t)_{∈ R}n_{, y(t)}

∈ Rℓ_{, and u(t)}

∈ Rm_{. Or, in discrete time,}

x(k + 1) = Adx(k) + Bdu(k)

y(k) = Cdx(k) + Ddu(k) (1.2)

The discrete-time model is the main interest of this thesis, since modern digital computers work with sampled data. In the sequel, we shall drop the subscripts c, d; since whether a model is in discrete or continuous time can be directly dis-tinguished from whether a differential or a difference state equation is used. Note that first principle based modeling often results in nonlinear models in continu-ous time. These can in fact be simplified to discrete time LTI models as in (1.2), following two steps:

• linearizing the nonlinear continuous time model at an equilibrium point; • discretizing the continuous time linear model at a given sampling rate. In reality, systems are usually subject to external disturbances, process and measurement noise, and faults. The model of (1.2) shall therefore be modified to incorporate these “extra inputs”.

x(k + 1) = Ax(k) + Bu(k) + Ef (k) + F w(k)

y(k) = Cx(k) + Du(k) + Gf (k) + v(k) . (1.3) Here w(k) ∈ Rnw _{and v(k)} _{∈ R}nv _{respectively represent the process and}

mea-surement noise. f (k) ∈ Rnf _{represents the fault inputs. A, B, C, D, E, F, G are}

bounded real matrices with appropriate dimensions.

Instead of incorporating the fault inputs f (k) in (1.3), another popular way in modeling actuator and sensor faults is modifying the B, C, D matrices (Gao and Antsaklis 1991,1992;Noura et al. 2000;Richter et al. 2007,2008). Total actuator faults can be modeled by annihilating the corresponding columns in B, i.e. for the total failure of the j-th actuator:

Bj ← 0 × Bj, where Bj∈ Rn×1. (1.4)

Total sensor faults can be modeled by annihilating the corresponding row in C matrix, i.e. for the i-th sensor:

Ci← 0 × Ci, , where Ci ∈ R1×n. (1.5)

Partial actuator/sensor faults are modeled via scaling the corresponding column/row of B/C matrix by a factor 0 < α < 1, i.e.

Bj ← α × Bj, (1.6)

Ci ← α × Ci. (1.7)

In fact, these models are special cases of the generic model of (1.3). To see this, note that the partial fault of the j-th actuator defined by (1.6) is equivalent to setting E = (α− 1)Bjand f (k) = uj(k) in (1.3), where uj(k) represents the control signal

(17)

to the j-th actuator. Due to this equivalence, we shall only consider the model (1.3) in this thesis, whose generality is further verified in the following detailed example.

1.2.2 Modeling faults in a three-tank system

h

H

Figure 1.3:Nominal configuration of the three-tank system

Under nominal condition, only the left (T1) and the middle tank (T2) are used.

The right Tank T3 and Pump P2 are not in use and act as redundant hardware.

The main aim of the two tanks used is to provide a continuous water flow q2to

a consumer. The water level in the middle T2has therefore to be maintained at a

level h2 ∈ [9, 11]cm. The reservoir tank T1is filled by Pump P1up to a nominal

water level of h1 = 50 cm. Water flows between the tanks can be controlled by

several valves_{V12H, V12L, V23H, V23L}. All valves can only be completely opened

or completely closed. A continuous positioning is not possible. For the nominal case, Valves_{V12L, V23H, V23L} are closed and not in use. Valve V12H is used to

control the water level in tank T2.

(18)

the Toricelli law, (1.8).

qij = Cij· sgn(hi− hj)·

q

|hi− hj| (1.8)

where Cijis a constant depending on the geometry of the connecting pipe and the

valve; hi, hj are the water levels above the connecting pipe. The change of water

volume V in a tank can be calculated as ˙

V = A_{· ˙h =}Xqin−

X

qout, (1.9)

wherePqinis the sum over all water inflows andPqoutis the sum over all water

outflows of the tank. A is the cross-section area of the cylindric tank and h the water level in the tank.

The details in modeling the whole system are given in Blanke et al. (2003);

Heiming and Lunze (1999). The complete model is nonlinear and time variant in the sense that the flow directions are determined by the level differences in the water tanks. We shall only focus on a linearized model around a nominal equilibrium point. We shall also include the leakage fault in T1to this model. This

linear model illustrates how actuator and (or) sensor faults can be incorporated into the linear state-space model in the form of (1.3). Note that in the nominal case, only the valve V12H is opened. Tank 3 is disconnected from the other two

tanks. The water level in Tank T1and T2is characterized by h1 ≥ hH ≥ h2. The

water hence flows from T1to T2, leaks out of T1, and supplies to the customer

from T2. In this case, one derives the mass balance equations as follows.

˙ h1 = 1 A(qP1 − q12H− qL), (1.10) ˙ h2 = 1 A(q12H− q2), (1.11) where qP1 = CP1uP1, q12H = C12H √ h1− hH, qL = CL√h1, and q2 = C2√h2. An

equilibrium is characterized by the constant water level in both tanks; i.e. ˙h1 =

˙

h2= 0. Denote an equilibrium by the superscript∗. Then at an equilibrium point,

one has q∗ P1 = q ∗ 12H+ q∗L ⇒ u∗P1 = C2 C_P1 p h∗ 2+CC_P1L p h∗ 1, q∗ 12H = q2∗ ⇒ h∗2 = C2 12H C2 2 (h ∗ 1− hH). (1.12) Denote α , C12H2 2AC2√h∗2 , β , CL 2A√h∗ 1 , ξ , C2 2A√h∗ 2 , and λ , CP1 A . The linearized

model thus takes the following form

d dt h1 h2 = −α − β 0 α −ξ h1 h2 + λ 0 uP1 = −α 0 α −ξ h1 h2 + λ 0 uP1+ 1 0 (_−βh1). (1.13)

(19)

Note that the leakage actually offsets the state matrix by −β 0 0 0 , which is therefore a plant fault. However, one can treat this as an additive fault to the model in a manner as in the second line of (1.13). This is in fact the case when one does not have the prior knowledge on the physical model of the leakage.

Denote the sampling time by Ts, the continuous time model, (1.13), will be

discretized as follows (Astr ¨om and Wittenmark 1984˚ ), h1(k + 1) h2(k + 1) = ζ1 0 φζ1+ ψζ2 ζ2 | {z } Ad h1(k) h2(k) + λ " 1−ζ1 α φ α(1− ζ1) + ψ ξ(1− ζ2) # | {z } Bd uP1(k) + Bd(− β λh1(k)) | {z } fd(k) , (1.14)

where ζ1 = e−αTs, ζ2 = e−ξTs, φ = _ξ−αα , and ψ = _α−ξα . (1.14) is obviously in the

form of (1.3).

Another fault in the three-tank system as defined inBlanke et al.(2003); Heim-ing and Lunze(1999) is the blockage of Valve V12H, when it is closed; i.e. no water

is supplied to the supply tank T2in the nominal configuration (q12H= 0). With the

same line of derivation, it is easy to show that the discrete-time linearized model around the equilibrium point as defined before can be written as,

h1(k + 1) h2(k + 1) = 0 0 0 ζ2 | {z } Ad h1(k) h2(k) + λ " 1−ζ1 α φ α(1− ζ1) + ψ ξ(1− ζ2) # | {z } Bd uP1(k) = ζ1 0 φζ1+ ψζ2 ζ2 | {z } Ad h1(k) h2(k) + λ " 1−ζ1 α φ α(1− ζ1) + ψ ξ(1− ζ2) # | {z } Bd uP1(k) + " 1−ζ1 α φ α(1− ζ1) +ψ−1ξ (1− ζ2) # | {z } Ed (αh1(k)) | {z } fd(k) , (1.15) with the same ζ1, ζ2, φ, ψ as given in (1.14). (1.15) also fits into the framework of

(1.3).

1.2.3 LPV state-space models with faults

The drawback of linearizing a nonlinear model at a certain working point is that the linearized model becomes invalid when the working point changes. One way to avoid linearization is to keep the nonlinear dynamic equations, but describe them in a special form; e.g. in a linear parameter varying structure (Wu 1995). Un-like an LTI model, time dependent scheduling parameters are present in an LPV

(20)

model (Shamma and Athans 1991;Wu et al. 1996;Wu 1995). The LPV description preserves the linear structure, but the state space coefficients are time-dependent. Specifically, we shall consider the following discrete time LPV model with ad-ditive faults under the affine form

x(k + 1) =

nµ

X

i=1

µ(i)_k A(i)x(k) + B(i)u(k) + E(i)f (k) + K(i)e(k) (1.16)

y(k) = Cx(k) + Du(k) + Gf (k) + e(k) (1.17)

where x(k) _{∈ R}n_{, y(k)}

∈ Rℓ_{, and u(k)}

∈ Rm_{. f (k)}

∈ Rnf _{represents the}

ad-ditive fault signals. The innovation signal e(k) is a zero-mean white noise with a non-singular covariance of EET_{. µ}(i)

k ∈ R are the scheduling parameters (or

local weights) of the local models, {A(i)_{, B}(i)_{, E}(i)_{, K}(i)

}. nµ denotes the

num-ber of these local models. A(i), B(i)_{, C, D, E}(i)_{, G, K}(i)_{are bounded real matrices}

with appropriate dimensions. The affinity, i.e. the linear combination of the ap-propriate local models times the scheduling signals, characterizes the parameter variation of the LPV model. The scheduling vector is supposed to be measurable at every sample k. We hereby claim that µk∈ Pµ, wherePµis a compact polytope.

The time varying state-space matrices can now be defined as,

A(µk) , Ak = nµ

X

i=1

µ(i)_k A(i). (1.18)

We shall now show a physical example of modeling nonlinear dynamics in such an LPV structure.

1.2.4 Modeling a DC motor with an unbalanced disc

A DC motor can be accurately described by a linear model (Pillay and Krishnan 1989). However, it becomes nonlinear, if one mounts an additional mass to the rotation disc, as depicted in Figure1.4. The mass distribution of the homogenous disc becomes in-homogenous. The mathematical description of the system can be decomposed to an electrical and a mechanical part. Table1.1contains the value and the nomenclature of the parameters.

On the first hand, using Kirchhoff’s Voltage Law,

Lm· ˙I(t) = u(t) − Ki· ω(t) − Rm· I(t), (1.19)

where I(t)[mA] is the current, u(t)[V ] the control input voltage and ω(t)[rad/sec] is the rotation speed of the disc. The latter variable connects the electric part to the mechanical; i.e.

J _{· ˙ω(t) = K}i· I(t) − b · ω(t) + M · g · l · sin(θ(t)), (1.20)

(21)

! L _R Ki U M.g.sin( ) M.g.cos( ) J, m l i

Figure 1.4:The schematic view of a DC motor with an unbalanced disk

Table 1.1:Nominal parameters of the DC motor

Parameter Value

Motor torque constant Ki= 53.6· 10−3N m/A

Motor resistance Rm= 9.50Ω

Motor impedance Lm= 0.84· 10−3H

Complete disk inertia J = 2.2· 10−4_{N m}2

Friction coefficient b = 6.6_{· 10}−5_{N ms/rad}

Additional mass M = 0.07 kg

Mass distance from the disk center l = 0.042 m

where the shaft angle is θ[rad]. Exact linear parametrically varying model can be derived if one replaces the nonlinear term with an appropriately chosen and mea-sured parameter. Let the scheduling parameter be defined by µ(t) = sin(θ(t))/θ(t). Since the angle is measured, the scheduling parameter can be evaluated. Besides, since_{| sin(θ(t))/θ(t)| ≤ 1, µ(t) is contained in a convex set, denoted as P}µ.

Combining the two parts into a single dynamics by measuring the angle and its derivative, one gets the following continuous time exact LPV state space form, where the state vector is x(t) =θ(t) ω(t) I(t)T.

 ˙ω(t)˙θ(t) ˙ I(t)   = "₀ ₁ ₀ 0 −b/J Ki/J 0 _−Ki/Lm −Rm/Lm # + " ₁ ₀ ₀ M gl/J 0 0 0 0 0 # µ(t) ! "_θ(t) ω(t) I(t) # + " ₀ 0 1/Lm # ·u(t), (1.22)

(22)

y(t) = 1 0 0 0 1 0  ω(t)θ(t) I(t)   . (1.23)

With the fast electrical dynamics neglected, the system becomes second order and affine in µ(t) with unchanged output equation (x(t) =θ(t) ω(t)T),

˙_θ(t) ˙ω(t) = 0 1 0 _−1/τ + 0 0 M gl/J 0 µ(t) θ(t) ω(t) + 0 Km/τ u(t), where τ = RmJ

Rmb+K2i is the modified time constant of the linear part; Km=

Ki

Rmb+Ki2.

The LPV model is full state observable and controllable. This can easily be checked by deriving the rank of the appropriate observability and controllability matrices or the invariant subspaces (Bokor and Balas 2004). Moreover, the pair (C, A(µ(t))) is quadratically detectable/stabilizable with parameter varying ob-server/controller feedback and simple Lyapunov solution matrix by solving the appropriate Linear Matrix Inequality (LMI) conditions on the vertices of the pa-rameter polytope Pµ (Wu 1995). The latter condition is of high importance to

describe the system in a stable, parameter varying and closed loop observer form. Rewriting Eqs. (1.22-1.23) and (1.24) in a continuous time and compact form:

˙x(t) = A(µ(t))x(t) + Bu(t),

y(t) = Cx(t). (1.24)

InToth et al.(2007) among the alternative ways of creating a discrete-time equiv-alence of (1.24), the rectangular method has been selected mostly because it pre-serves the nature of the scheduling parameter µ and does not increase the system dimension; i.e.

x(k + 1) = Ad(µk)x(k) + Bdu(k),

y(k) = Cx(k), (1.25)

with the sampling time Ts. Here, Ad(µk) = In+ Ts· A(µk), Bd = Ts· B and In

is the state dimensional identity matrix. Consequently, the model with additive fault inputs can be modeled as

x(k + 1) = Ad(µk)x(k) + Bdu(k) + Edf (k),

y(k) = Cx(k) + Gf (k). (1.26)

Respectively for additive actuator and sensor faults, Ed = Bd, G = 0 and Ed =

0, G = Iℓ.

1.3 Model based fault detection and identification

1.3.1 Approaches for LTI systems

(23)

u(k) Residual

Generator

Plant Fault

Estimation f (k)ˆ

y(k) r(k)

Figure 1.5:The block diagram of an FDI filter, where u, y, r, ˆf represent respec-tively inputs and outputs of the plant, residuals, and estimated faults.

detection can be regarded as a residual generation problem. A residual can be defined in Definition1.2.

Definition 1.2 (Residual (Isermann and Ball´e 1997)) A fault indicator, based on a deviation between measurements and model-equation based computations.

According to the survey of Frank(1990), there are four different approaches in residual generation.

• the parity space approach (PSA), • the dedicated observer approach, • the fault detection filter approach, • the parameter identification approach.

PSA is a model-based statistical detection method monitoring a sequence of residuals along a receding horizon. Model based PSA for LTI systems has been intensively studied by many authors (Chow and Willsky 1984;Desai and Ray 1981;Ding et al. 1999;Frank 1990;Gertler and Singer 1990;Patton et al. 1989). See Section3.1in this thesis. The robustness of PSA against modeling errors and/or disturbances has been addressed by Chen and Patton(1999);Frank(1994);Lou et al.(1986). In a recent work (Zhang and Ding 2007), PSA has also been extended to periodic systems. It is also worth mentioning here that the moving horizon esti-mation (MHE) approach, which is perhaps more popular in the predictive control literature (Muske et al. 1993;Rao et al. 2001), has also been applied in FDI based on the receding horizon principle (Tyler et al. 2000).

Similar to MHE, observer based and fault detection filter design approaches are related to state estimation (Frank 1990). In observer based approaches, Luen-berger observers or Kalman filters are designed based on a nominal plant model. In fault detection filter approaches, the filter is optimized for special purposes, e.g. for isolating a single fault. In fact, the state observer based design has been one of the main stream research directions in the FDI literature since 1970’s, according to the survey ofFrank(1990). This approach yields a residual observer based on the state-space model of a plant. Due to the recent progress in_H∞control and linear

matrix inequality (LMI) techniques (Boyd et al. 1994;Zhou et al. 1996), the latest interest along this track has been focused on improving the robustness of an FDI filter against model uncertainties in anH∞ framework (Frisk and Nielsen 2006;

(24)

Mangoubi 1998), and improving the sensitivity of such a filter in the presence of perturbations in a mixed_H₋/_H_∞framework (Ding 2008;Henry and Zolghadri 2005;Wang et al. 2007;Wei and Verhaegen 2009).

The approaches reviewed above mainly deal with additive faults. As far as multiplicative faults are concerned, the parameter identification approaches can be opted, since this category of faults change the parameters of a system ( Gustafs-son 2001). This approach relies on recursive system identification (Ljung 1987). Another recent development along this track lends itself to the literature of mul-tiple model approaches, originally developed for radar detection (Li and Jilkov 2005). This approach models a plant under various conditions (nominal or faulty) in a (convex) model set (Hallouzi et al. 2008;Li et al. 2005), and detects and iso-lates faults according to the magnitudes of the online estimated model weights (Hallouzi et al. 2006;Zhang and Li 1998).

Fault identification concerns locating the fault occurred in a system. This is done by comparing the non-zero residual to pre-computed fault vectors, and the closest one in a statistical sense indicates the location of the fault (Blanke et al. 2003;Gustafsson 2001,2007). Faults can also be isolated, if their magnitudes can be estimated. In the existing literature, fault estimation has been postulated as a filtering problem (Stoustrup and Niemann 2002), as an MHE problem (Qin and Li 2001;Tyler et al. 2000), and as a system inversion problem (Hou and Patton 1998;

Sundaram and Hadjicostis 2008;Szigeti et al. 2001).

1.3.2 Approaches for LPV systems

For general nonlinear systems, developing an FDI approach becomes more in-volved. Taking the parity space approach as an example. Three difficulties appear due to nonlinearity (Bokor 2009). First, the parity relation becomes nonconstant and nonlinear. Second, there is no general realization theory, which can lead to a general parity space approach for nonlinear systems. Besides, for general nonlin-ear systems, even if state variables can be eliminated, a residual becomes a non-linear function of not only faults, but also I/Os and their derivatives (Bokor 2009). This nonlinearity denies a closed-form expression of the statistics of the residual. A linearization of the nonlinear residual generator around the nominal param-eters is therefore needed (Zhang et al. 1998). The classical statistical detection theories (Basseville and Nikiforov 1993) can then be applied to the neighborhood of the nominal parameters. This solution for general nonlinear systems is how-ever prone to the errors in the linearization step (Zhang et al. 1998). In view of the difficulties associated with the general nonlinearity, special model structures have to be considered to simplify the problem.

Recently, the FDI for LPV systems is gaining more popularity in the literature (Balas et al. 2002;Bokor and Balas 2004;Bokor 2009;Kulcs´ar 2006;Wei and Ver-haegen 2008). The LPV description preserves the linear structure, but the state space coefficients are time-dependent. This brings an advantage over a general nonlinear model description; i.e. the residual generated by the input-output rep-resentation of an LPV model has a linear dependence on stochastic disturbances.

(25)

No linearization is required in this case, provided the scheduling parameters are fault free. We shall postpone the detailed discussion about the statistical detection approach of LPV systems to Section3.4of this thesis. Generally, the FDI scheme for an LPV system can be represented by the diagram in Fig.1.6, where the resid-ual generator and fault estimator become time dependent due to the parameter varying nature of the LPV model.

u(k) Residual Generator LPV Plant Fault Estimation f (k)ˆ y(k) r(k) µ(k)

Figure 1.6:The block diagram of an FDI filter for LPV systems, where the LPV model is defined in Section1.2.3.

1.4 Model based control reconfiguration

A control system that can tolerate faults is designed based on two philosophies: • designing a controller, which is robust to any faults that can be foreseen to

occur in the system;

• designing an “adaptive” controller, which adjusts itself based on online FDI information.

These two types of FTCs are also respectively referred to as passive and active designs in the literature (Blanke et al. 2003).

Passive FTCs are usually postulated as robust control synthesis problems in theH∞ framework. The parameter changes due to multiplicative faults can be

treated as model parameter uncertainties, and hence naturally falls into the do-main of robust optimal control (Liao et al. 2002;Niemann and Stoustrup 2002;

Zhou et al. 1996;Zhou and Ren 2001). On the other hand, additive faults can be treated as disturbances, whose effects on the controlled outputs can be rejected in anH∞sense (Campos-Delgado et al. 2005) or in the sense of Lyapunov

stabil-ity (Mhaskar et al. 2006). Recently,Kanev et al.(2003,2004);Kanev(2004) have extended the passive design based on the classical robust control to possibly non-convex parameter sets, using the technique of randomized algorithms (Calafiore and Dabbene 2006;Tempo and Dabbene 2001).

However, in reality, potential faults in a system can hardly be foreseen. The passive designs may cover as many faulty situations as possible, by enlarging either the uncertain parameter set for multiplicative faults or the norm of the un-known additive faults. The price to pay is the increased conservativeness and reduced closed-loop performance. This is a drawback of the passive approaches,

(26)

and motivates the research on the active approaches. Active FTCs rely on online FDI information, instead of pre-assumed faulty scenarios.

One direction in the active design is the model following reconfigurable con-trol. The approaches ofGao and Antsaklis(1991,1992) aim at finding a state feed-back gain such that the state matrix of a closed-loop faulty system is as close as possible to that of the nominal system, in terms of a Frobenius norm. The state-space model of a faulty system is assumed to be known. Recently, similar ap-proaches have been developed to specifically handle actuator and sensor faults (Richter et al. 2007,2008). This approach aims at recovering the I/O behavior of a plant under actuator faults, through matching the Markov parameters of the nominal and faulty plant. The reconfigured closed-loop plant “hides” the actua-tor faults, and therefore referred to as the “fault-hiding” property. The solution for LTI systems has recently been extended to piecewise affine systems (Richter et al. 2008). These methods require that the actuator or sensor faults are correctly de-tected and isolated. Similar to the model following or the fault hiding techniques, eigenvalue assignment methods have also been reported in the literature (Zhang and Jiang 1999,2000). These methods try to achieve the same eigenvalues of the closed-loop faulty system as those of the nominal one, through state-feedback.

The methods reviewed above all require the knowledge of the faulty model. If the faults to tolerate come from actuators or sensors, the faulty model can be online determined by isolating the fault, and annihilating the corresponding columns in the input matrix (see Eq. (1.4)) or rows in the output matrix (see Eq. (1.5)). As long as component (or multiplicative) faults are concerned, the faulty model cannot be determined, unless the plant parameters are identified recur-sively online. This then motivates another direction in the active FTCs, i.e. the adaptive control approaches (Goodwin and Sin 1984; Astr ¨om and Wittenmark˚ 1994). InTao et al. (2002a,b), model following reconfigurable control methods have been addressed. Unlike the approach of Gao and Antsaklis (1991,1992);

Richter et al.(2007,2008), the parameters of the closed-loop transfer function are identified online. Parameter adaptation based reconfigurable schemes are also developed byBoˇskovi´c and Mehra(2005);Ichtev et al.(2002);Wang et al.(1997);

Zhang et al.(2004);Zhang and Jiang(2001).

Despite the numerous publications in FDI and CR, an integrated design still re-mains a hard topic. The challenges are attributed to the imperfect FDI information contaminated with delays and errors (Mahmoud et al. 2003). Recent contributions have been reported byBoˇskovi´c and Mehra(2005);Noura et al.(2000);Zhang et al.

(2004);Zhang and Jiang(2001).

1.5 Classical designs based on identified models

The FDI and CR approaches reviewed in the previous sections are based on mod-els. In this section, we take a closer look into the procedures in the classical model-based control designs. Model-model-based FDI filters can be designed through similar procedures, and shall be omitted for brevity.

(27)

1.5.1 Building a model by system identification

The state space models (1.1) and (1.2) can be built by first principles. However, such a modeling is often a time consuming and expensive task (Hjalmarsson 2009). This has stimulated the development of system identification techniques (Ljung 1987;S ¨oderstr ¨om and Stoica 1988;Verhaegen and Verdult 2007). In this sec-tion, we shall limit the modeling discussions to identification based approaches.

As long as identification is concerned, there exist two different approaches, prediction error methods (PEM) (Ljung 1987;S ¨oderstr ¨om and Stoica 1988) and subspace model identification (SMI) (Verhaegen and Verdult 2007).

Prediction error methods aim at identifying the following parameterized mod-els,

y(k) = G(z, θ)u(k) + H(z, θ)e(k). (1.27) Here y, u, e are respectively outputs, inputs, and white noise. z is the forward shift operator. θ represents the set of parameters to identify. G and H are transfer functions, according to the chosen model structure, e.g. ARX, ARMAX, output error, Box-Jenkins, and etc. The design parameters in a prediction error method are therefore the structure and order of G and H. For instance, an ARMAX model takes the following structure,

y(k) = b1z−1+· · · + bnbz−n b 1 + a1z−1+· · · + anaz−na u(k) + 1 + c1z−1+· · · + cncz−n c 1 + a1z−1+· · · + anaz−na e(k). (1.28)

In this case, na, nb, nc need to be chosen. The estimation of the coefficients is a

nonlinear optimization problem. However, for multivariable systems, one en-counters the difficulty of parameterizing multioutput models, because the inter-nal structure of these models become much richer (Ljung 1987). For instance, in a mulivariable ARMAX model, G, H change to matrix fraction descriptions, e.g.

G_{(z, θ) = A}−1_(z)B(z), H_{(z, θ) = A}−1_(z)C(z),

whose elements are transfer functions. In prediction error methods, the identifica-tion of G, H can be carried out channel by channel, followed by a model reducidentifica-tion step (Patwardhan and Shah 2005). To reduce the complexity of such an approach, connection is established with subspace identification (Ljung 1987).

Subspace methods tackle the difficulties in identifying multivariable systems by estimating either the range space of the extended observability matrix ( Verhae-gen and Verdult 2007) or the unknown sequence of its states (Chiuso 2007a;van Overschee and de Moor 1994) of the following state space model,

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

y(k) = Cx(k) + Du(k) + e(k). (1.29)

Depending on the parameterization of (A, B, C, D, K), this model corresponds to different structures of (1.28) (Verhaegen and Verdult 2007). It is hence the identifi-cation data, instead of a user, that decides the model structure fitting it best. Since

(28)

the methodologies to be developed in this thesis are targeted for multivariable systems, we shall only consider subspace identification methods.

Recently, a new predictor-based subspace identification (PBSID) approach has been developed byChiuso(2007a,b) for consistent identification from data mea-sured in closed-loop plants. Due to this consistency, we will especially consider the PBSID approach in this thesis. For simplicity, we shall take D = 0 in (1.29) to introduce this method in rest of this subsection.

Since the innovation e(k) is unknown, it can be replaced by y(k)_{− C ˆx(k). A} closed-loop observer thus results from (1.28); i.e.

ˆ x(k + 1) = (A_{− KC)} | {z } Φ ˆ x(k) + Bu(k) + Ky(k), y(k) = C ˆx(k) + e(k). (1.30)

The following data equation is then a consequence of applying (1.30) in predicting the outputs along a horizon of L sampling instants, where i = 0,_{· · · , N − 1 and} B =B K. This equation can be regarded as a lifted model of the state-state form (1.30) with L I/O pairs. Note that there are only two tuning levers in this data equation; i.e. the past and future horizon s, L.

     y(k + i) y(k + i + 1) .. . y(k + i + L_{− 1)}      | {z } ˆ y[k+i,k+i+L) =      C CΦ .. . CΦL−1      | {z } OL ˆ x(k + i)+      0 CB 0 .. . ... CΦL−2_B _CΦL−3_B _{· · · 0}      | {z } TL,z        u(k + i) y(k + i) .. . u(k + i + L− 1) y(k + i + L− 1)        | {z } z[k+i,k+i+L) +      e(k + i) e(k + i + 1) .. . e(k + i + L− 1)      | {z } e[k+i,k+i+L) . (1.31) Let the following data matrices be defined.

Yk+N +L−2 k = ˆ y[k,k+L) yˆ[k+1,k+L+1) · · · ˆy[k+N −1,k+N+L−1) , Ek+N +L−2_k = e[k,k+L) e[k+1,k+L+1) · · · e[k+N −1,k+N+L−1) , Xk,N = x(k)ˆ x(k + 1)ˆ · · · ˆx(k + N − 1), Zkk+N +L−2 =       

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

..

. ... ...

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

       = z_[k,k+L) z_[k+1,k+L+1) _{· · · z}_{[k+N −1,k+N+L−1)} .

(29)

k_{− s}

past

k_{− 1 k}

future

k + L_{− 1}

Figure 1.7:Illustration of the “past” and the “future” horizon in SMIs.

Eq. (1.31) can now be extended to the following matrix-form data equation.

Ykk+N +L−2=OL· Xk,N+TL,z· Zkk+N +L−2+ Ek+N +L−2k . (1.32)

Since the state sequence Xk,N is usually unknown, (1.32) actually formulates

a nonlinear optimization problem. Based on (1.30), the productOL· Xk,N is

con-structed from the measurable I/O data to tackle this difficulty; i.e. OL· Xk,N=    CΦs .. . CΦs+L−1    · Xk−s,N | {z } bx +      CΦs−1_B _CΦs−2_B _{· · ·} _CB CΦs_B _CΦs−1_B _{· · ·} _CΦB .. . ... ... CΦs+L−2_B _CΦs+L−3_B _{· · · CΦ}L−1_B      | {z } Hs,z · ·       

u(k− s) u(k − s + 1) · · · u(k − s + N − 1) y(k_{− s) y(k − s + 1) · · · y(k − s + N − 1)}

..

. ... ...

u(k_{− 1)} u(k) _{· · ·} u(k + N_{− 2)} y(k_{− 1)} y(k) _{· · ·} y(k + N_{− 2)}

       | {z } Zk+N −2 k−s . (1.33) By assuming Φ stable (Ljung 1987), the bias bx can be neglected by choosing a

large s. Combining respectively (1.31) and (1.33) and (1.32) and (1.33), one derives the following data matrices,

ˆ

y[k,k+L) ≈ Hs,zz[k−s,k)+TL,z· z[k,k+L)+ e[k,k+L), (1.34)

Ykk+N +L−2 ≈ Hs,zZ_k−sk+N −2+TL,z· Zkk+N +L−2+ Ek+N +L−2k . (1.35)

Here the “≈” sign is due to ignoring bx. An extra time window of s sampling

instants is introduced. Note that the first column of_Z_k−sk+N −2and_Z_kk+N +L−2is re-spectively z_[k−s,k)and z[k,k+L). These two windows are hence respectively called

the “past” and the “future” horizon, as illustrated in Fig. 1.7. The value of s and L can be respectively chosen as s_{≫ 1, L ≥ 1. Eq. (}1.35) now becomes a linear LS problem to estimate the unknown parametric matrices_Hs,z,TL,z.

As a summary, the following steps are required by the PBSID approach. Step1. Solve the linear LS problem (1.35) forHs,z,TL,z.

(30)

Step2. Estimate_OL· Xk,N by (1.33), with bxignored. Compute the following

sin-gular value decomposition (SVD),

Hs,z· Z_k−sk+N −2=U1 U2 Σ1 Σ2 · VT 1 VT 2 . (1.36)

Here Σ1 contains the n dominant singular values of Hs,zZ. Up to users’

choice, the number n is the selected model order. The state sequence can now be estimated by ˆXk,N = Σ1V1T. Then by solving a least-squares

prob-lem with ˆXk,Nand the measured I/Os, the matrices A, B, C, K in (1.30) can

be identified, subject to a similarity transformation.

We shall show the effectiveness of this approach by an example.

Example 1.1 (Subspace identification) For simplicity, consider identifying the state space model of the following SISO plant,

(1− 0.8z−1_{)y(k) = (z}−1_{− 0.9z}−2_{)u(k) + e(k).} _(1.37)

Here e(k) represents the measurement noise of a zero mean and a variance of 10−4_{. Let}

us be given a sequence of 1000 I/O samples measured from the closed-loop plant of (1.37)

controlled by an initial controller u(k) = r(k)−0.5·y(k), where r(k) is a white noise with

a zero mean and a variance of 0.01. The horizons were chosen as s = L = 6. The model order indeed turned out to be 2, according to the singular values ofHs,zZ. The identified

poles are depicted in Fig. 1.8. The result verifies that the more general state-space model structure of (1.29) identified in the SMI holds for an ARX model structure, even though

no ARX model is explicitly parameterized before being fitted to data.

−1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Real axis Imaginary axis

Figure 1.8:Identified poles of the plant (1.37) by subspace methods. Cross: real; “x”: identified.

(31)

1.5.2 Designing a controller based on an identified model

With a state space model, e.g. (1.29), a wide variety of control schemes can be opted to design an optimal control strategy; e.g. LQG, H∞, and MPC. In this

thesis, we are particularly interested in predictive control, because it is connected to the data equation (1.31) defined in the subspace identification, and can hence be integrated into the unified framework of identification, control, fault detection, and experiment design to be developed in this thesis.

In fact, as MPCs have become standard industrial practice, there are numerous software toolboxes available to do the synthesis and implementation. With the help of these tools, the main efforts in applying an MPC do not lie in formulating the predictive control law itself any more, but rather in modeling the dynamic process.

MPCs (Camacho and Bordons 2004;Maciejowski 2000) regulate future outputs to referenced trajectories. The following future Np-step ahead output predictor is

hence required to be built from a state space model, e.g. (1.29). ˆ y[k,k+N p)=     C CA .. . CANp−1      | {z } ONp x(k) +      D CB D .. . ... CANp−2_B _CANp−3_B _{· · · D}      | {z } TNp      u(k) u(k + 1) .. . u(k + Np− 1)      | {z } u_[k,k+Np) . (1.38) The decision variables are u(k),_{· · · , u(k + N}c− 1); i.e. setting u(k + Nc) =· · · =

u(k + Np− 1) = u(k + Nc− 1), where Nc≤ Npis the control horizon. See Fig.1.9

for an illustration of the horizons.

k-1 k+ 1 k+ Nc k+ Np y(k+ i)^ u(k+ i) r(k+ i) prediction horizon control horizon k

(32)

The block elements in the Toeplitz matrix_TNp are the so-called Markov

pa-rameters of the plant. Denote the sequence of Markov papa-rameters in the sequel as M_u→y ₌_CANp−2_B _{· · · CB D}_.

The initial state, x(k), can be estimated by a Kalman filter; i.e. ˆ

x(k + 1) = (A− KC)ˆx(k) + (B − KD)u(k) + Ky(k),

y(k) = C ˆx(k) + Du(k) + e(k). (1.39)

This filter is in the same form as the closed-loop observer form of (1.30). The identified state space matrices by an SMI can be used as an approximate Kalman filter (van Overschee and de Moor 1994;Verhaegen and Verdult 2007).

MPCs optimize the cost function,JMu→y,(1.29)paramererized by Mu→yand the

state-space model (1.29), subject to I/O constraints.

minu[k,k+Nc) JMu→y,(1.29)(u[k,k+Nc), ˆx(k))

s.t. u[k,k+Nc)∈ U

ˆ

y[k,k+Np)∈ Y

. (1.40)

Here_{U, Y are convex I/O sets. The objective function takes the following form,} JMu→y,(1.29) =kr[k,k+Np)− ˆy[k,k+Np)k

2

Q+ku[k,k+Nc)k

2

R. (1.41)

Here Q, R_{0. r}[k,k+Np)represent the future reference signals. Out of the optimal

u[k,k+Nc), only u(k) is applied to the plant at the current time instant k, which is

known as the receding horizon principle.

The design and implementation procedures of the classical MPC approach is illustrated in Fig.1.10.

Assume L = Np. Note that Eq. (1.33) indicates

ONpx(k)ˆ ≈ Hs,z· z[k−s,k). (1.42)

Rewriting (1.38) in terms of the Kalman filter form (1.39) and substituting (1.42) into the resulted equation, one can see that the output predictor exactly resembles the data equation (1.35). Here, ˆx(k) can be regarded as the Kalman filter esti-mate, since the identified K matrix is an approximate steady-state Kalman gain. Based on this fact,Favoreel et al.(1998,1999) have suggested to use the identified Hs,z,TL,zat the first step of an SMI1 to formulate the output predictor (1.35). It

can be seen now that the common data equation (1.35) represents a striking resemblance

in identification and control analysis.

Because of their connection with subspace identification, the control laws de-veloped byFavoreel et al.(1998,1999) are called subspace predictive control, or

1_{Note that in the work of}_{Favoreel et al.}₍₁₉₉₈_,₁₉₉₉_{), the Markov parameters in H}

s,z,TL,zare the

(33)

x(k + 1) = Ax(k) + Bu(k) + Ke(k) y(k) = Cx(k) + Du(k) + e(k)

M_u→y₌_CANp−1_B _{· · · CB D} minu_[k,k+Nc) JMu→y,(1.29) s.t. u[k,k+Nc)∈ U ˆ y[k,k+Np)∈ Y measure identify compute MPC control past I/Os

Figure 1.10:The concept of classical MPCs. The solid and dashed lines respec-tively represent the offline and online procedures.

SPC. Similar to (1.41), an SPC boils down to the following optimization problem, minu[k,k+Nc) JHs,z,TL,z(u[k,k+Nc), z[k−s,k))

s.t. u[k,k+Nc)∈ U

ˆ

y[k,k+Np)∈ Y

. (1.43)

Here the objective function_J_Hs,z,TL,ztakes the same form as in (1.41); except that

ˆ

y[k,k+Np) is parameterized byHs,z,TL,z, and depends on the past I/Os z[k−s,k),

instead of the Kalman filter estimate ˆx(k). The SPC scheme can be depicted as in Fig.1.11. Note that the term “SPC” should be distinguished from “statistical process control”, which also appears in fault detection literature (Kourti and Mac-Gregor 1996).

Note that in the existing SPC approaches, the unknown Markov parameters are simply replaced by their estimates, without taking the estimation errors into account. We shall show through the following example that when the errors are large, this “certainty equivalence” solution fails.

Example 1.2 (Subspace predictive control) Consider again the SISO plant of (1.37)

given in Example1.1. Increase the variance of e(k) to 0.16; while keeping the other pa-rameters and experiment conditions unchanged. We shall illustrate that the parameter estimates contain large errors, in terms of the identified magnitude response of the plant as shown in Fig.1.12.

We consider the closed-loop SPC approach ofFavoreel et al.(1999), with s = L =

Np = 6. Set Q = I6, R = 0.001I6, where I6 is a 6× 6 identity matrix. The control

objective is to regulate the plant output to 0, from the initial value of 100. The result is shown in Fig.1.13. Clearly, the SPC could not stabilize the closed-loop plant.

(34)

Hs,z,TL,z minu_[k,k+Nc) JHs,z,TL,z s.t. u[k,k+Nc)∈ U ˆ y[k,k+Np)∈ Y measure identify SPC control past I/Os

Figure 1.11:The concept of SPCs. The solid and dashed lines respectively repre-sent the offline and online procedures.

10−3 10−2 10−1 100 0.5 1 1.5 2 Frequency (Hz) Magnitude

Figure 1.12:Identified magnitude response of the plant (1.37) from noisy I/O measurements. Solid: real; Dashed: identified.

This example brings a clear message that the SPCs have to be robustified against parameter estimation errors. In the literature of adaptive control, such a robust so-lution is known as cautious control (Campi and Prandini 2003;Ohrn et al. 1995¨ ;

Zhang and Guo 2006). We shall introduce this design concept with the approach developed byZhang and Guo(2006).2

Let the prediction horizon of (1.35) be one. The following one-step ahead

out-2_{Note that the approach of}_{Zhang and Guo}₍₂₀₀₆_{) is developed for SISO plants and with only one}

step ahead output prediction. For MIMO plants and multiple steps of output prediction, the analysis becomes more involved, because of the matrix form of the Markov parameters and the propagation of their errors along the prediction horizon. See Section2.4of this thesis for details.

(35)

5 10 15 20 25 30 35 −200 −150 −100 −50 0 50 100 150 200 samples

Figure 1.13:Outputs of the closed-loop plant (dashed curve) under the control of the SPC approach of Favoreel et al.(1999), when the parameter estimates contain a large error.

put predictor results, ˆ y(k + 1) =CΦs−1_B _CΦs−1_K _{· · · CΦB CΦK CK} | {z } ξs ·φ[k−s,k+1)+ CBu(k), where φ[k−s,k+1) = uT(k− s), yT(k− s), · · · , uT(k− 1), yT(k− 1), yT(k) T . A “certainty equivalence” design, minimizing (ˆy(k + 1)_{− r(k + 1))}2_{, is given by}

u(k) =₋ξˆsφ[k−s,k+1)− r(k + 1)

CB . (1.44)

Here r(k + 1) represents the output reference; CB and ˆξsare respectively the

iden-tified CB and ξs. A cautious solution to the problem of minu(k)E(ˆy(k + 1)− r(k +

1))2takes the following form (E is the expectation operator),

u′(k) =₋ CBξˆsφ[k−s,k+1)− r(k + 1) + P2s 2s−1φ[k−s,k+1) (CB)2_{+ P} CB . (1.45)

Here PCB represents the second last diagonal element of the covariance matrix P

of the identified parameters. P_2s−12s stands for the second last row of P .3_{Now, we} shall illustrate this cautious design by the following example.

Example 1.3 (Cautious control) Consider again the SISO plant of (1.37) given in

Ex-ample1.1, where the variance of e(k) is set to 0.16, as in Example1.2. ξsand CB were

identified with the past horizon set to s = 6. Both the certainty equivalence control (1.44)

and the cautious control (1.45) were applied, to regulate the plant output to 0, from the

initial value of 100. The 2-norm of the covariance matrix P was found to be 0.03. PCB

turned out to be 0.0183. The I/O signals are shown in Fig.1.14.

Clearly, the magnitude of the cautious control signal was always smaller than that of 3_{See Section}_2.4_{for the detailed discussion about the error covariance matrix.}

(36)

the certainty equivalence control signal, after both the algorithms started from the 6-th sampling instant. This is due to the extra penalty put on the cautious control due to the variance of the estimated parameters. This can be seen from (1.45), whose denominator

yields q (CB)2_{+ P} CB > CB. 0 50 100 150 200 −100 0 100 200 output 0 50 100 150 200 −100 −80 −60 −40 −20 0 input samples

Figure 1.14:Outputs of the closed-loop plant under the control of both the “cer-tainty equivalence” control (1.44) (dash-dotted curve) and the cau-tious control (1.45) (dashed curve).

Note also from Example1.3that, although the cautious control has robustness against parameter estimation errors, it suppresses the magnitude of the control inputs. This is a bad news for adaptive control, where the parameter convergence depends on the persistence of excitation of the inputs (Campi and Prandini 2003). On the other hand, in terms of predictive control, the cautious design based on Markov parameters has an advantage over that based on state-space matrices. To see this, let us consider the uncertainties in Φ and B and assume no uncertainty in C, K; i.e. Φ = Φ0+ ∆Φ, B = B0+ ∆B. In this case,

CB = C(B0+ ∆B) = CB0+ C∆B,

CΦB = C(Φ0+ ∆Φ)(B0+ ∆B) = CΦ0B0+ CΦ0∆B + C∆ΦB0+ C∆Φ∆B

| {z }

∆(CΦB)

.

∆(CΦB) is quadratic in ∆Φ, ∆B, and does not necessarily define a convex set. Similarly, ∆(CΦτB), τ > 1 are of even higher orders in ∆Φ, ∆B, and hence non-linear and nonconvex. Therefore, the output predictor based on these state-space matrices becomes (highly) nonlinear in ∆Φ, ∆B, and hence leads to a nonlinear and nonconvex predictive control law.

As a summary, the model-based and subspace predictive control schemes are compared in Table1.2.

(37)

MPC SPC

Parameters state-space matrices Markov parameters

needed (A, B, C, D, K) (D, CΦi_{B K}_{, i}

≥ 0)

Subspace two LS problems, SVD, one LS problem

Identification order selection

Uncertainty in nonlinear, linear,

output predictor nonconvex convex

Advantages control of hidden states, few design parameters, stability/performance analysis online reconfiguration Table 1.2:Comparisons between MPC and SPC.

1.5.3 Fault detection based on an identified model

With a state space model, e.g. (1.29), a wide variety of fault detection schemes can be opted, as reviewed in Section1.3. In this thesis, we are particularly interested in parity space approach, because of its residual generation in a moving horizon, where an output predictor can also be connected to the data equation (1.31) as explained below.

Classical PSA considers the state space model (1.3), which leads to the follow-ing lifted input-output relation in a horizon of Npsamples,

y[k,k+N_ p)=     C CA .. . CANp−1      | {z } ONp x(k) +      D CB D .. . ... CANp−2_B _CANp−3_B _{· · · D}      | {z } TNp      u(k) u(k + 1) .. . u(k + Np− 1)      | {z } u_[k,k+Np) +      G CE G .. . ... CANp−2_E _CANp−3_E _{· · · G}      | {z } TNp,f      f (k) f (k + 1) .. . f (k + Np− 1)      | {z } f_[k,k+Np) +      0 CF 0 .. . ... CANp−2_F _CANp−3_F _{· · · 0}      | {z } TNp,w      w(k) w(k + 1) .. . w(k + Np− 1)      | {z } w_[k,k+Np) +      v(k) v(k + 1) .. . v(k + Np− 1)      | {z } v_[k,k+Np)

(38)

Or in a compact form,

y[k,k+Np) = ONpx(k) +TNpu[k,k+Np)+TNp,ff[k,k+Np)

+ TNp,ww[k,k+Np)+ v[k,k+Np)

(1.46)

The following residual is then defined rp,k = P· y[k,k+Np)− TNpu[k,k+Np) = P_·_ONpx(k) +TNp,ff[k,k+Np)+TNp,ww[k,k+Np)+ v[k,k+Np) , (1.47) where P defines a projection, which yields P· ONp = 0, P· TNp,w = 0 and P·

TNp,f 6= 0. Although eliminating the influence of the unknown initial state x(k)

and unknown disturbance w[k,k+Nc), the projection reduces the degree of freedom

(DOF) in the residual rp,k, and may hence hide possibly useful information therein

(T ¨ornqvist 2006).

Assume the disturbances w(k), v(k) are white noise. The state space form (1.3) has the following Kalman filter

ˆ

x(k + 1) = (A− KC)ˆx(k) + (B − KD)u(k) + (E − KG)f(k) + Ky(k) y(k) = C ˆx(k) + Du(k) + Gf (k) + e(k).

The statistics of the innovation e(k) are determined by those of w(k), v(k) (Ljung 1987). See also Section2.2for the details. The PSA residual generator (1.47) can hence be rewritten in terms of this Kalman filter form. Similar to the SPC scheme introduced in the last subsection, the unknown productONpx(k) can be approxi-ˆ

mated byHs,zz_[k−s,k). The elimination ofONpˆx(k) by projecting rp,konto the left

null space ofONp can hence be avoided. This thesis tries to use this idea to

de-velop new fault detection and identification algorithms. Due to their connection with the data equation (1.35) in the SMI, these new algorithms are called “Fault detection and Identification approaches Connected to Subspace Identification”, or FICSI.

1.6 Scope of the thesis

Through the discussions in the last section, it is clear now that subspace identifi-cation, predictive control, and fault detection can be unified, in the sense that they can all be formulated based on the output predictor, i.e. the data equation in the form of Eq. (1.35) and the approximation of_ONpˆx(k) byHs,z· z[k−s,k) as in Eq.

(1.42). This unified framework can hence be illustrated by Fig.1.15. The parame-ters needed in the output predictor can be identified from data by solving a linear LS problem. This indicates that this unified framework links to the physical world via data, instead of first principles, and is hence entitled with “data driven”. Data driven designs are recently gaining more popularity, due to the latest progress in sensor and data storage technologies (Bars et al. 2008). According toQin(2009), data driven methods are relatively easier to implement in real processes of rather

(39)

Output Predictor Identification Predictive Control Fault Detection

Figure 1.15:The unified framework of identification, predictive control, and fault detection. The intersection indicates dependence.

large scale, when compared with other methods based on rigorous process mod-els.

Obviously, the quality of the data determines the convergence of parameter estimation, and hence the performance of data driven designs. Here, the data quality does not refer to the signal to noise ratio (SNR), but to the information that they contain on the dynamics of the system. “Informative” data have to be gener-ated by input signals, which can persistently excite (PE) the system (Ljung 1987;

S ¨oderstr ¨om and Stoica 1988). Guaranteeing the PE condition is the major concern of the classical input design approaches. However, the classical designs can only be performed offline and for a single input channel. In this thesis, we aim at de-veloping an online approach for MIMO systems, such that the reconfiguration of the SPC scheme can be carried out online. Moreover, PE input design is especially important for cautious adaptive control, where the control inputs are suppressed in the presence of parameter estimation errors.

As a summary, the scope of this thesis is to develop data driven control and filter designs, which maximize as much as possible the information that can be retrieved from the data without the user to interfere. This design objective is real-ized for LTI systems by the unified framework of identification, control, FDI, and persistently exciting input design, where the parameter identification is a single linear least-squares problem. The thesis also aims at developing the unified iden-tification, predictive control, and FDI design to LPV systems. The overall scheme of the data driven designs developed in this thesis is illustrated in Fig.1.16.

1.7 Organization of the thesis

Chapter2of this thesis provides the solutions of the data driven control design, i.e. the closed-loop subspace predictive control (SPC). The concept of the SPC is

Data Driven Fault Tolerant Control: A Subspace Approach