Predictive control: A lecture course given in the Aerospace Engineering Faculty TU Delft

Iseries 03:

Control and Simulation 05

Predictive Control

A Lecture Course Given in the Aerospace

Engineering Faculty TU Delft

J. M.

Maciejowski

Predictive Control

A Lecture Course Given in the

Aerospace Engineering Faculty TU Delft

Bibliotheek TU Delft

111111111111

C 3021888

2392

337

4

Series 03: Control and Simulation 05

. ;

..

.

.

Predictive Control

A Lecture Course Given in the

Aerospace Engineering Faculty TU Delft

J.M. Maciejowski

Published and distributed by:

Delft University Press Mekelweg 4

2628

CD Delft The Netherlands Telephone

+31 (0)152783254

Fax

+31 (0)152781661

e-mail: DUP@DUP.TUDelft.NL by order of:

Faculty of Aerospace Engineering Delft University of Technology Kluyverweg 1

P.O.

Box

5058

2600 GB

Delft The Netherlands Telephone

+31 (0)152781455

Fax

+31 (0)152781822

e-mail: Secretariaat@LR.TUDelft.NL website: http://www.lr.tudelft.nl/

Cover: Aerospace Design Studio,

66.5

x

45.5

cm, by:

Fer Hakkaart, Dullenbakkersteeg

3, 2312

HP Leiden, The Netherlands Tel.

+ 31 (0)71 51267 25

90-407-1714-1

Copyright ©

1 998

by Faculty of Aerospace Engineering

All rights reserved.

_.1111 ..... + . ... .....

No part of the material protected by this copyright notice may be reproduced or utilized in any farm or by any means, electronic or

mechanical, including photocopying, recording or by any information storage and retrieval system, without written permission from the publisher: Delft University Press.

Preface

These notes were written for a course of 14 two-hour lectures on Model Predictive Control given to students at the Technical University of Delft in November and December 1997.

I would like to thank Bob Mulder and Hans van der Vaart for inviting me to spend a sabbatical term at the Faculty of Aerospace Engineering in Delft, for making all the arrangements, and for making my stay at Delft both pleasant and interesting. I would also like to thank Ton van den Boom for giving me feedback on the first few chapters, and Rob de Vries for some valuable help with chapter 4.

Contents

Preface

1 Introduction to Predictive Control 1.1 Motivation . . . . 1.2 General Optimal Control Formulation 1.3 The 'Receding Horizon' Idea . . . . 1.4 Terminal Constraints Ensure Stability 1.5 Problems . . . ..

2 A Standard Formulation of Predictive Control 2.1 General Features of Constrained Predictive Control . 2.2 A Standard Formulation

2.3 Prediction . . . . . . . .

2.3.1 No disturbances, FUll state measurement . 2.3.2 Constant output disturbance

2.3.3 Using an observer 2.4 Problems .. . . . .

3 Solving Predictive Control Problems 3.1 The MPC Tooibox . . . . 1 1 2 4 8 12 13 13 15

20

20

23 24 28 31 31

CONTENTS

3.2 Unconstrained Problems . . . . .. 34 3.2.1 Measured State, No Disturbances . 34 3.2.2 Formulation as a Least-Squares Problem . 36 3.2.3 Structure of the Unconstrained Controller 38

3.2.4 Estimated State 40

3.3 Constrained Problems 42

3.3.1 Formulation as a QP Problem . 42

3.3.2 Controller Structure 44

3.4 Problems .. . . 48

4 Step Response and Transfer Function Formulations 51

4.1 Step and Pulse Response Models 51

4.1.1 Step and Pulse Responses 51

4.1.2 Relations to State-Space Models 54 4.1.3 Prediction Using Step Response Model. 56 4.1.4 State-Space Models From Step Responses 58

4.2 Transfer Function Models 65

4.2.1 The Basics

.

.

..

65

4.2.2 Prediction Using Transfer Functions 68 4.2.3 Prediction with a Disturbance Model. 71

4.2.4 The GPC Model .. . . 77

4.2.5 State-Space Interpretation . 78

4.2.6 MultivariabIe Systems 81

4.3 Problems

...

.

.

..

.

..

.

82

CONTENTS

5.1 What Are We Trying To Do?

5.2 Some Special Cases . . . . . . 5.2.1 Effect of Control Weighting 5.2.2 Mean-Level Control 5.2.3 Dead-Beat Control 5.2.4 'Perfect' Control

5.3 Frequency-Response Analysis

5.4 Disturbance Models and Observer Dynamics . 5.4.1 Disturbance models

.

.

.

5.4.2 An implementation issue. 5.4.3 Observer dynamics 5.5 Set-point Pre-filtering 5.6 Problems

...

....

6 Stability 6.1 Terminal Constraints . 6.2 Infinite Horizons .. .

6.2.1 Infinite horizons give stability .

6.2.2 Feasibility . . . . . . . . . . . .

. . . . . . . .

.

.

....

6.2.3 Constraints and infinite horizons - stable plant 6.2.4 Constraints and infinite horizons - unstable plant 6.3 Fake Algebraic Riccati Equations

6.4 Using the Youla Parametrization 6.5 Constraints and Feasibility 6.6 Problems . . . . . . . . . . 85 89 89 89 90 91 94 96 96 101 103 104 108 111 111 112 112 114 114 116 118 122 124 126

CONTENTS

7 Robust Predictive Control 127

7.1 Formulations of Robust Control . 127

7.1.1 Norm-bounded uncertainty 127

7.1.2 Polytope uncertainty . 128

7.2 LQG /LTR Approach . 128

7.3 LMI Approach 129

7.3.1 Overview 129

7.3.2 Robustness without constraints 130

7.3.3 Robustness with constraints . 134

7.4 Problems

.

....

...

..

..

..

.

137

8 Perspectives 139

8.1 Parametrising the Control Input 139

8.2 Zones, Funnels and Coincidence Points . 139

8.3 Spare Degrees of Freedom . 140

8.3.1 Ideal resting values . 140

8.3.2 Multiobjective formulations 141

8.3.3 Fault tolerance . 142

8.4 Constraint Management 144

8.5 Nonlinear Internal Models 146

8.6 Non-quadratic Costs 146

8.7 Concluding Remarks 146

List of Figures

1.1 Receding Horizon Predictive Control. Parts band d show the situ-ation at time t, while parts a and c show it at time t+TI, with TI = 1. (Taken from the book by Soeterboek. [46]) .. . . .. . . .. . . 7

2.1 Paper machine with headbox 18

2.2 Output disturbance. 23

2.3 A State Observer .. 25

3.1 Structure of controller with no constraints and full state measurement. 39 3.2 Structure of controller with no constraints and state observer. 41 3.3 Structure of controller with constraints and state observer.. . 46

4.1 Step response of 3-input, 2-output distillation column. Solid curves: original responses. Dotted curves: responses of 8-state model. 63 4.2 The first 15 singular values of 11.250 for the example. . . . . . 64 4.3 Spectral density of crew disturbance, modelled by (4.108) with p = 0.98. 72 4.4 Generation of v(klk) - one interpretation. 74

5.1 Two degree of freedom feedback system. . . 87 5.2 Mean-level control: set-point change on rl. 93 5.3 Mean-level control: set-point change on r3. 94 5.4 Default (DMC) observer: response to step disturbance .. 95 5.5 Deadbeat observer: response to step disturbance . .. . . 96

LIST OF FIGURES 5.6 Deadbeat set-point response. Step demand on rl .. . . . 97 5.7 'Perfect' controller: Set-point response. Step demand on rl. 98

5.8 'Perfect' controller: Set-point response. Step demand on r3. 99 5.9 Singular values of S(z) (solid line) and T(z) (broken line) for

mean-level control of the headbox. . . . . 101

5.10 A 'proportional + integral' controller. 102

5.11 Predictive controller, unconstrained case. 104 5.12 Predictive controller, simplified. . . . . 105

5.13 Predictive controller, simplified again. 106

6.1 Youla parametrization of all stabie feedback systems (if the plant is stabie). . . . . .. . . .. . . 123

Chapter

1

Introduction to Predictive

Control

1.1

Motivation

Most courses on Control Engineering deal with contral of linear systems. But most systems to be controlled are nonlinear - even if the dynamics can be approximated by linear equations, they have actuators which saturate, or can move at limited rates, and their outputs cannot reach values outside limited ranges. Also, most courses on Control deal with Single-Input, Single-Output (SISO) systems. But most systems to be controlled have several inputs and outputs. Control design techniques for linear systems have been extended to multivariable systems [26, 52], but most of these do not say what to do about nonlinear phenomena - the designer has to deal with them in some unsystematic method.

Predictive Control is a method of controlling systems which deals with nonlinearities which arise from constraints, and also deals naturally with multivariable systems. It has become very widely and successfully used in certain sectors of the process industries, particularly the petrochemical sector. In these industries, maximum profit is usually obtained from a process by running it as close as possible to some constraints: pressures or flow rates as high as possible, for example. Predictive control operates by performing optimization subject to constraints, and as aresult allows operators to run their processes much closer to constraints than would be possible with conventional linear controllers (such as 'proportional and integral' controllers, for example). Increasing efficiency or profitablity by 1% or 2% can be worth millions of dollars per year on a single plant in the petrochemicals sector, and remarkably short pay-back periods have been reported for predictive control installations.

Another use for Predictive Control is as a basis for Adaptive Control. Since pre-1

2 CHAPTER 1. INTRODUCTION TO PREDICTIVE CONTROL dictive control uses an explicit 'internal model' of the plant being controlled, it is natural to think of updating the model if the plant changes, and in this way obtain adaptive contro!. This mode of using adaptive control has been particularly popular with academic researchers, and it has also been applied to some extent in industry. As we shall see, Predictive Control involves solving optimization problems on-line. It does not appear suitable, therefore, for aerospace applications. There is the problem of computation time, which can be considerable (but not a problem in many process industries, where control update times of several minutes are acceptable), and there is the problem of certification, since it is difficult to say exactly what a predictive controller is going to do in a given situation. The problem of computation time will gradually disappear as computer speeds continue to increase. But even now, there are aerospace applications for which today's computers would be fast enough to use predictive contro!. Large flexible spacecraft, for example, with very large solar panel arrays (such as the projected International Space Station) have very slow dynamics and do not require fast update rates. Even in aircraft flight control, it is possible to think of using predictive control for 'higher-level' functions such as flight management, in which the constraints due to safety and comfort are combined with air traffic control requirements to find trajectories which are 'efficient' in some sense [39].

Terminology:

Various names have been used for Predictive Contro!. 'Model Predictive Contro!' or 'MPC', 'Model Based Predictive Control' or 'MBPC', 'Constrained Predictive Control' , 'Receding Horizon Control' , 'Generalized Predictive Control' or 'GPC', all mean the same thing. There are also trade names of proprietary predictive control products, such as 'Dynamic Matrix Control' or 'DMC', and 'Model Algorithmic Control' or 'MAC'.

1.2

General Optimal Control Formulation

The idea of posing control problems as .problems of constrained optimization is not new. In fact all the 'classical' theory of Optimal Control, as developed between, say, 1955 and 1970, was diven by problems of constrained optimization arising out of the needs of the aerospace industry, particularly by military needs in flight and missile control, and by the problems of launching, guiding and landing space vehicles. In some ways this theory solved an extremely wide range of problems. Suppose that the 'plant' being controlled has an input vector u and a state vector x, and has nonlinear behaviour governed by the vector differential equation

dx

dt = !(x, u).

1.2. GENERAL OPTIMAL CONTROL FORMULATION 3

function') which has the form

V(x, u, t) =

foT

f(x(t), u(t), t)dt

+

F(x(T))

where f(x(t), u(t), t) is some function which is never negative, and that the control input is constrained to be in some set u(t) E U.

This is a very general formulation, and just about every sensible problem can be represented in this form by using suitable functions j, F and f.

Example 1.1 Suppose the plant is a rocket which is to be transjerred from the launch site - state Xo - to a low-earth orbit - state x I - using as little fuel as possible, and the vector of inputs inc/udes the juel flow rate as the first element, Ul. This problem can be represented by setting

f(x(t), u(t), t) = IU1(t)1 and

F(x(T))

=

{oo

if x(T)

#

xI

o

ij x(T)

=

xI

The vector junction f depends on the dynamics of the rocket. The optimization will be subject to the constraint 0 ::; Ul(t) ::; umax , where Umax is the maximum possible fuel flow rate.

In principle, the solution to such general problems is known: you have to solve the so-called Hamilton-Jacobi-Bellman equation [5J:

a

°

. (

a

o(

-a _t

v

(x,t) = mmH x,u'-a V x,t)) uEU x

where H(x, u, À) = f(x, u) +Àf(x, u), with the boundary condition VO(x, T) = F(x).

Once you have solved this, you choose the control input signal as the one which minimises H:

uO(x, t) = arg minH(x, u, /Ja Vo(x, t)). uEU X

Note that this is a feedback controllaw, because it depends, at time t, on x(t), which itself depends on what has happened over the whole interval from the initial time to the current time. (Note: You do not have to know where all this comes from. The point is to see that it is nasty.)

Unfortunately, it is virtually impossible to solve this partial differential equation in most cases! A few particular problems have been solved (of the type: 'Minimse fuel (or time) needed to get from A to B'), but that is all. 80 it is necessary to consider some more specific problems in order to get solutions in practice.

4 CHAPTER 1. INTRODUCTION TO PREDICTIVE CONTROL

One way to go is to assume that the plant is linear, so that the function

f

(x, u)

takes the special form:

:i;

=

A(t)x(t)

+

B(t)u(t)

and that the functions f. and F are quadratic functions:

f.(x(t), u(t), t)

F(x(T))

X(t)T Q(t)x(t)

+

u(t)T R(t)u(t) x(T)TSx(T)

where Q(t), R(t) and S are square, symmetric, positive-definite matrices. In that

case the Hamilton-Jacobi-Bellman equation simplifies to an ordinary differential

equation (the Riccati equation) which can be solved, and which leads to a linear

(but time-varying) feedback law:

u(t) = F(t)x(t)

where F(t) is a matrix (the 'state feedback' matrix) which depends on the solution of

the Riccati equation [5, 22]. Probably more than 95% of the development of optimal

control has been devoted to going in this direction. It leads to a very complete and

powerful theory, which has been applied very widely - particularly in its 'dual'

form, for the filtering problem, in which case it leads to the theory of Kalman filters. But it has nothing to say about constraints or other nonlinearities.

The other way to go is down the road of Predictive Control.

1.3 The 'Receding Horizon' Idea

The difficulty of the general optimal control problem arises because it is a function

optimization problem. Instead of ordinary calculus, it needs the Calculus of

Vari-ations, to find the optimal input function among all possible ones. (For a very nice

account of how this theory developed see [47].)

The main idea of Predictive Control is to avoid this difficulty by restricting the set of possible inputs to such an extent that the optimization is performed only over a

finite set of parameters or 'decision variables', rather than over a set of functions.

Most commonly this is done by cutting up time into discrete intervals, optimising

over a finite horizon (T), and hence over a finite number of intervals, and assuming

that the input signal is held constant during each interval.

Imagine we are at time t, and we have measurements from the plant, so we know where we are starting from. For the time being, let us assume that we can measure

the state x, so we know everything ab out the plant. Formulate an optimal control

problem over some future horizon {t + T : 0

:s

T

:s

T}:

1.3. THE 'RECEDING HORJZON' IDEA 5

Mini-Tutorial: Quadratic Forms

Expressions like x(t)TQ(t)x(t) and u(t)TR(t)u(t), where x(t), u(t) are vectors and Q(t), R(t) are symmetric matrices, are called Quadratic Forms. They are just compact repres-entations of certain quadratic functions in several variables. This is best illustrated by an example.

Example 1.2 Suppose that

x

[

:~

]

and Q

[

10 4

4 3 .

]

Then x(tfQ(t)x(t) _{[Xl X2}

₁

[

~O ~][:~]

10xr + 8X1X2

+

3x~.

Notice that we do not use any special notation for vectors and matrices. You will have to teil from the context whether a symbol represents a number, vector, or matrix. (But we will usually use upper-case letters for matrices.)

If xT Qx

>

0 for every x except x = 0, then this quadratic form is called positive-definite.

It can be shown that aquadratic form is positive definite if and only if all the eigenvalues of the matrix involved in it are positive, in which case we write Q

>

o.

The matrix Q in the example has both eigenvalues positive. So the quadratic form is positive-definite. You can try substituting any values for Xl and x2 (even negative ones), and you will see that the value is always positive. Note that this happens even if Xl = 0 or

X2 = 0, provided that they are not both zero.

If Q

=

I then xTQx

=

xT X

=

Ilx11

2, the squared 'length' of the vector, or squared Euclidean norm. Even if

Q =I-

I, it is possible to think of aquadratic form as the square of a 'weighted norm' (if Q

>

0), and the notation xTQx = Ilxll~ is of ten used.

The only other thing we will need to know ab out aquadratic form is how to find its gradient. Let V = xTQx, but now consider the general case where X has n elements, x = [Xl, X2, . .. , x n f, and Q is an n x n matrix. Then the gradient of V is given by:

\7V

["-,,,-, ... ,.,,--1

8V 8V 8V

UXl UX2 UX n

2xTQ

(We are not going to prove this here.) Note that we have defined the gradient to be a row vector. Sometimes it is defined to be a column vector, in which case it is given by 2Qx.

6 CHAPTER 1. INTRODUCTION TO PREDICTIVE CONTROL

where T is now fixed, and we have simplified things by assuming that f. depends only on the predicted state, x(t+rlt) and on the input û(t+r), which we think now (at time t) we will apply later (at time t + r). These predicted values will probably be different from the actual state x(t + r) and the actual input u(t + r) which will occur at time t + r. Note that we have simplified things by assuming that f. does not dep end explicitly on t.

Minimising V(x(t),u) over the set of admissible input signals {u(t + r) : 0 ::; r ::;

T, u(t

+

r) E U} gives us a proposed optimal input signal {u°(t

+ rit) :

0 ::; r ::; T}. The next key idea is that we only apply an initial portion of this input signal, say for a time interval TI, and then we throw away the rest of the optimal input signal {u(t

+rlt) :

TI

:S

r

:S

T}, measure the current state X(t+TI) again, and re-solve the optimization problem by minimising V (x( t

+

TI), u), getting a new optimal input signal {uO(t

+

TI

+ rit +

TI) : 0 ::; r ::; T}. This process is repeated indefinitely. Figure 1.1 shows the basic idea. This figure shows the set-up we are going to adopt in this course, namely dividing up time into regular intervals, and keeping the input signal constant during these intervals. Notice how the prediction horizon (Hp 'time steps' long in the figure) slides forward as time goes on, or 'recedes'. Hence the term Receding Horizon.

Figure 1.1 should be examined closely, because it contains a number of important features. Firstly, time has been discretised, as has already been said. Secondly, notice the distinction between 'real time', measured by t, and 'time now', measured by k, which slides along in real time. Thirdly, the control signal u changes only a finite number of times, and is assumed to remain constant af ter that. This is crucial in practice, because it means that the optimization problem we will have to solve will involve only afinite number of'decision variables' {û(k+ilk) : i = 0, 1, ... ,Hu -I}, with Hu::; Hp.

Although at each step an open-loop problem is solved, it is always based on the latest measurement of the state, so that overall a feedback control policy is implemented. So any unexpected disturbances which occur, or modelling errors, are taken account of.

Since we will work with discrete time, we move from a differential equation model of the plant to a difference equation model:

x(k

+

1) = f(x(k),u(k)) and we make a corresponding change to the cost function:

Hp

V(x(k), u) =

L

f.(x(k

+

ilk), û(k

+

i - 11k))

+

F(x(k

+

Hplk)) i=1

1.3. THE 'RECEDING HORIZON' IDEA a b I k I k-l k+l

k~lkk~l

k+Hp past - - - - future I i i I I I I I t I - time t-l t+l

c

d I k I k-l k+l I k I k-l k+1 7 u u past - - - - future I I I I I t I t-l t+l - t i m e

Figure 1.1: Receding Horizon Predictive Control. Parts band d show the situation at time

t,

while parts a and c show it at time

t

+

TI, with TI = 1. (Taken from the book by Soeterboek. [46])

(and usually we will not use the penalty on the predicted terminal state F(x(k

+

Hplk)). Note that the function

e

depends on the predicted state at time k

+

i and the input one step earlier, at time k

+

i - 1, since u(k

+

i) has no effect on x(k

+

i)

according to the model.

The receding horizon idea requires V (x (k ), u) to be optimized on-line in real time, with the model (1.3) acting as a constraint, and other constraints due to equipment

and safety limitations. (1.3) and (1.3) are still much too genera! for this to be feasible, and we will focus on some very specific choices in the rest of the course.

8 CHAPTER 1. INTRODUCTION TO PREDICTIVE CONTROL

1.4 Terminal Constraints Ensure Stability

Predictive Control, using the receding horizon idea, is a feedback control policy.

There is therefore a risk that the resulting closed loop might be unstable. Even

though the performance of the plant is being optimised over the prediction horizon, and even though the optimization keeps being repeated, each optimization 'does not

care' about what happens beyond the prediction horizon, and so could be putting

the plant into such a state that it will eventually be impossible to stabilise. This is particularly likely to occur when there are constraints on the possible control

input signais. It is easy to construct an example where this happens, even without

constraints:

Example 1.3 Consider a discrete-time system consisting of two delays in series: xl(k

+

1) u(k)

x2(k

+

1) xl(k) or, using vector-matrix notation:

x(k) =

[~ ~]

x(k)

+

[~]

u(k)

Suppose the prediction horizon is Hp = 1, and the cost to be minimised is V(x(k),u(k)) xl(k

+

11k)2

+

6X2(k

+

11k)2 +4Xl(k

+

1Ik)x2(k

+

11k)

=

x(k

+

11kf

[~ ~]

x(k

+

11k)

(which is positive-definite). The predictive control problem is to find u(k) which will minimise V(x(k),u(k)).

Using the model, the predictions are xl(k

+ 11k)

=

u(k) and x2(k

+ 1

1k)

=

xl(k).

Substituting these into V (x( k), u( k)) gives

V(x(k),u(k)) = u(k)2

+

6Xl(k)

+

4u(k)Xl(k).

Ta find the minimum of this, dijJerentiate and set the derivative to zero: 8V

8u(k) 2u(k)

+

4Xl(k)

o

ifu(k) = - 2Xl(k).

1.4. TERMINAL CONSTRAINTS ENSURE STABILITY 9 or

x(k +

1)

=

[~2 ~]

x(k).

Clearly this is unstable, since the closed-Ioop transition matrix has an eigenvalue at -2, outside the unit circle.

You might guess that this problem arises in the Exarnple because the prediction horizon is too short, so that the control is too 'short-sighted'. This is correct, and it turns out that stability can usually be ensured by making the prediction horizon long enough, or even infinite. We will see that later in the course.

Another way of ensuring stability is to have any length of horizon, but to add a 'terminal constraint' which forces the state to take a particular value at the end of the prediction horizon. And it is surprisingly easy to prove this, even in a very general case, by using a Lyapunov function. For the purpose of proving stability it is enough to consider the case when the state is to be driven to the origin, from some initial condition.

Theorem 1.1 Suppose that Predictive Control is obtained for the plant x(k

+

1) = f(x(k), u(k))

by minimising the cost function Hp

V(x(k), u) =

L

l(x(k

+

ilk), û(k

+

i - 11k)) i=1

where l(x, u) 2 0 and l(x, u)

=

0 only if x

=

0 and u

= 0

, and l is 'decrescent',

subject to the terminal constraint

The minimisation is over the input signals {û(k

+

ilk) : i = 0,1, . .. , Hu - I}, with Hu

=

Hp (for simplicity), subject to the constraints û(k+ilk) E U and x(k+ilk) EX, where U and X are some sets. We assume that x

=

0, u

=

0 is an equilibrium condition for the plant: 0 = f(O,O). The receding horizon method is applied, with only uO(klk) being used from the optimal solution {uO(k

+

ilk) : i = 0, ... ,Hu -I}. Then the equilibrium point x

=

0, u

=

0 is stable, providing that the optimization problem is feasible and is solved at each step.

Proof 1.1 Let VO(t)

=

V(x(t), ua) be the optimal value of V which corresponds to the optimal input signal uO, as evaluated at time t. Clearly VO(t) 20, and VO(t)

=

0 only if x(t)

= 0 -

because if x(t)

= 0 then the optimal thing

to do is to set u(Hi) = 0

10 CHAPTER 1. INTRODUCTION TO PREDICTNE CONTROL for each i. We are going to show that V°(t

+

1) ~ VO(t), and hence that VO(t) is a Lyapunov function for the closed-Ioop system.

As usu al in stability proofs, we will assume that the plant model is perfect, 50 that the predicted and real state trajectories coincide: x(t+i)

=

:î:(t+ilt) ifu(t+i)

=

û(t+ilt).

With this assumption we have Hp

VO(t

+

1) = mJn

L

l(x(t

+ 1 +

i), u(t

+

i)) i=l

~

{~I(X(t

+ ;),

u(t - 1

+ ;)) -

I(x(t +1),.(t))+ l(x(t

+

1

+

Hp), u(t

+

Hp))}

~ -l(x(t

+

1), uO(t))

+

VO(t) + min {l(x(t + 1 + Hp), u(t

+

Hp))}

u

since the optimum is certainly no worse than keeping the optimal solution found at time t, which will take us up to time t

+

Hp, then doing the best we can for the final step. But we have assumed that the constraint x(t

+

Hp) =

°

is satisfied, since the optimization problem was assumed to be feasible, 50 we can make u(t

+

Hp)

=

°

and stay at x = 0, which gives

min{l(x(t+1+Hp),u(t+Hp))} = 0.

u

Since l(x(t),uO(t)) ~ 0, we can conclude that VO(t

+ 1)

~ VO(t). So VO(t) is a Lyapunov function, and we conclude by Lyapunov's theorem that the equilibrium x

=

0, u

=

°

is stabIe.

This sounds too good to be true. It seems too easy to guarantee stability. The problem is with the assumption that the optimization problem has a solution at eaeh step, and that it ean be solved at eaeh step. General eonstrained optimization problems ean be extremely diffieult to solve, and just adding a terminal constraint may not be feasible. But we will use the idea of this proof several times later, in more realist ie situations.

Example 1.4 If we look again at Example 1.3 and try to get stability by adding the terminal constraint x(k

+

1)

=

0, we get an infeasible problem, because that could only be aehieved if x1(k) = 0, which will not be troe in general. (In general you need at least 2 steps to drive a 2-state linear system to the origin.) Problem 5 shows that for this example stability is obtained by increasing the prediction horizon to Hp = 2, even without a terminal constraint.

1.4. TERMINAL CONSTRAINTS ENSURE STABILITY 11

Mini-Tutorial: Stability and Lyapunov Functions

Equilibrium: The equation x(k + 1) = f(x(k), u(k)) has an equilibrium (or 'fixed point')

at state Xo and input Uo if Xo = f(xo, uo).

Note that we can always introduce a change of coordinates z(k) = x(k) - xo, v(k)

=

u(k) - uo, so that the equilibrium is at (0,0) in the new coordinates:

z(k + 1)

=

x(k + 1) - Xo

o

f(z(k) + xo,v(k) + uo) - f(xo,uo)

g(z(k), v(k)) g(O,O)

So we can always assume that an equilibrium is at (0,0).

Stability: For nonlinear systems one has to consider the stability of a particular equilibrium

point, rather than the system. (Some equilibria may be stable, others unstable.)

An equilibrium point (0,0) is stabIe (in the sense of Lyapunov) if a small perturbation of

the state or input results in a 'continuous' perturbation of the subsequent state and input

trajectories. More precisely, given any € > 0, there is a ó > 0 (which depends on € > 0 and

the system) such that if II[x(Of,u(0)T111

<

€ then II[x(kf,u(k)TlIl

< ó for all

k

>

O.

It is asymptotically stabIe ifin addition II[x(k)T,u(k)Tlll-t 0 as k -t 00.

For closed-loop systems u(k) itself depends on x(k), so we could have written x(k + 1)

=

f(x(k)), where f now represents the closed-loop state equation.

Lyapunov's Theorem: If there is a function V(x, u) which is positive-definite, namely

such that V(x, u) 2 0 and V(x, u) = 0 only if (x, u) = (0,0), and has the ('decrescent')

property that

and if, along any trajectory of the system x(k

+

1)

=

f(x(k), u(k)) in some neighbourhood

of (0,0) the property

V(x(t + 1),u(t + 1))

:s:

V(x(t),u(t))

holds, then (0,0) is a stabie equilibrium point.

If, in addition, V(x(t),u(t)) -t 0 as t -t 00 then it is asymptotically stable.

12 CHAPTER 1. INTRODUCTION TO PREDICTIVE CONTROL

1.5

Problems

The first 3 problems can be skipped by those familiar with quadratic forms, Y'V, etc.

1. Use Matlab to find the eigenvalues of the matrix Q which appears in Example 1.2. Check that they are both positive.

2. (a) Suppose V(x)

=

9xi + 25x~

+

16x~. Find Q such that V(x)

=

xTQx. (b) Suppose V(x, u)

=

(5xi

+

2x~

+

x~)

+

(100ui

+

4u~). Find Q and R such

that V(x) = xTQx

+

uTRu.

(c) Are the functions V which appear in this problem positive-definite? 3. Check that the formula (1.1) for the gradient is correct for Example 1.2, by

working out the partial derivatives in the elementary way, and comparing the results.

4. Simulate the model and controllaw from Example 1.3 using Matlab, and check that the closed loop is unstable.

Recommended procedure: Create the state-space matrices a, b, c ,d corres-ponding to the closed-loop system, then create a Matlab object representing this closed-Ioop system using:

c1sys=ss(a,b,c,d,l);

(Use help ss to see what this does.) Then get the response from initial con

-dition xO= [1; 0] - or some other initial condition - by typing initial(c1sys,xO);

Note that these recommendations are suitable for use with Matlab version 5 only.

5. Consider the plant given in Example 1.3 again. But now suppose that the prediction horizon is Hp = 2, with the cost being

V(x(k), u(k))

=

x(k

+

1jkf

[~ ~]

x(k

+

Ijk)

+

x(k

+

2jk)T

[~ ~]

x(k

+

2jk) (a) Keep Hu = I, so that only u(k) is to be optimised, with the assumption

that u(k

+ 1)

= u(k). Show that the predictive controllaw is uO(k)

=

-~xl(k), and hence that the closed loop is stabie.

(b) Now let Hu = 2, so that u(k) and u(k+1) have to be optimised. By setting both derivatives (8V/8u(k)) and (8V/8u(k

+ 1)) to zero (or

Y'uV = 0 if you are happy with that notation) show that the predictive control law is uO(k)

= -~xl(k),

and hence that the closed loop is stabie.

Chapter 2

A Standard Formulation of

Predictive Control

2.1

General Features of Constrained Predictive Control

Since predictive control problems usually include constraints, the resulting control law is usually nonlinear. It is easy to understand why this is so. Suppose we have a gas-turbine engine with a 'hard' state constraint of the form xi(t)

<

Xi - such as 'Turbine temperature lower than 1000 °C'. If the turbine is running at 995°C and a disturbance d occurs which has the effect of moving xi away from Xi - the cooling oil flow rate increases for some reason, say - then the control system will react in a fairly 'relaxed' way, and will coordinate bringing xi back to its set-point with maintaining other objectives - correct turbine speed, power output, fuel economy,

etc. Now if a disturbance -d occurs instead, then by definition a linear controller would react in a similarly relaxed way, replacing a control signal u(t) by -u(t).

In contrast, a controller which is aware of the constraint will react in 'panic mode'. It 'knows' that there is no point in worrying about other objectives, becauseif the constraint is exceeded catastrophic failure is likely to result. It therefore reacts in a very different way to the disturbance -d than it would to disturbance +d - in our example, perhaps it will cut the fuel flow rate sharply and change the inlet vane angle to get maximum air flow. In a multi-engine aircraft it might even shut the engine down. Such behaviour is of course nonlinear.

We will see later that 'standard' predictive controllers in fact behave linearly so long as the plant is operating safely away from constraints, but nonlinearly when constraints are approached too closely. On the other hand, a standard predictive controller would not go so far as to shut down an engine, and it does not have the authority to do so. That kind of decision is usually handled by a separate higher-level 'supervisor', which may be a human operator or pilot.

14CHAPTER 2. A STANDARD FORMULATION OF PREDICTIVE CONTROL Although constrained predictive controllers are nonlinear, they are usually time-invariant. This means that there is a function h, such that the control signal can be expressed as u

=

h(x), namely it depends only on the current state, and not explicitlyon time (that is, it is not necessary to write u = h(x, t)). Of course this function 'exists' only in a mathematical sensej it would be impossibly difficult to write it down in 'closed-form' (as a formula). In order to see what u is, you have to feed x into an optimization algorithm and see what comes out. This time-invariance occurs so long as there is nothing in the problem formulation that depends explicitly on time: the model of the plant, the cost function and the constraints must all be independent of time.

We have to be rather careful about what this means exactly. First of all, our model

x(k

+

1)

=

f(x(k), u(k)) (2.1)

does not dep end on k except through x and u. Secondly in the cost function Hp

V(x(k), u)

=

L

e(x(k

+

ilk), û(k

+

i -11k), i) + F(x(k

+

Hplk)) (2.2) i=l

we can have explicit dependence of

e

on i but not on k. 80 you can attach a different weight to an error predicted to be 3 steps in the future than to an error predicted 1 step ahead. But when you solve the problem at the next k value, you must have the same pattern of weights. In a sense, i is a kind of 'dummy variabie' which only appears inside the control al!}}rithm, but is not seen from outside - in fact you can see this, because the sum Ei;l e( ... ) does not dep end on i. As time moves on, so k increases to k

+

1, but exactly the same problem i solved as was solved at time k.

80 the cost function does not dep end on time.

Thirdly, we have to be careful ab out constraints. Again we can have constraints depending on i but not on k. 80 IXj(k+ilk)1

<

Xj(i) is ok, but IXj(k+ilk)l

<

Xj(k) is not. More generally,

is ok, whereas x(k

+

ilk) E X(i) û(k

+

ilk) E U(i) x(k

+

ilk) E X(k) û(k

+

ilk) E U(k) (2.3) (2.4) (2.5) (2.6) is not. This can be arestriction if we want a time-invariant controller. For instance, if you are landing an aircraft, you may want constraints to become tighter as the aircraft approaches the runway. If you assume a certain speed profile, you can 'schedule' the constraint tightening as a function of time. Then you will have a time-varying predictive controllaw. However it may make more sense to schedule the constraint tightening as a function of distance from the runway, and make this

2.2. A STANDARD FORMULATION 15 one of the states. In this case the control law will be time-invariant (though not necessarily simpier!).

There is no necessary advantage to having a time-invariant controllaw. But we will see that having one may make analysis of its behaviour easier in some cases.

2.2

A Standard Formulation

Even though we have already simplified the general formulation of an optimal control problem to the constrained predictive control problem defined by (2.1)-(2.4), it is still too general and difficult to solve in real time. Further specialisation of the formulation is required.

From now on, we shall assume that the plant model is linear, that the cost function is quadratic, and that the constraints are in the form of linear inequalities. We shall also assume that that everything is time-invariant. Furthermore, we shall assume that the cost function does not penalise particular values of the input vector u(k), but only changes of the input vector, .ó.u(k), which are defined by

.ó.u(k) = u(k) - u(k - 1) (2.7)

To make the formulation useful in the real world, we shall drop the assumption that all the state variables can be measured, and that we only care about what the state variables are doing. Finally, up to now we have assumed that we always want to drive the plant to x = 0; it is more realistic to assume that we want some linear combination ofthe states, say z(k) = Czx(k), to follow a 'reference' trajectory

{rek)}. These assumptions coincide (almost) with those made in the great majority of the predictive controlliterature.

Let us formulate all this more precisely. We take the real input to the plant at time k to be u(k), the real plant state to be x(k), and the real measured output to be y(k). We use x(k

+

ilk) to denote our estimate, made at time k, of what the state will be at time k

+

i (and use similar notation for other variables ).

Plant Model:

States: x(k + 11k) = Ax(klk) + Bu(k) Measured outputs: y(klk)

=

Cyx(klk)

+

Dyu(k) Controlled outputs: i(klk) = Czx(klk)

+

Dzu(k)

(2.8)

(2.9) (2.10) We will usually assume that the real plant is governed by the same equations as the model, although this is not really true.

Cost function:

Hp Hu-l

V(x,u,r)=

L

lIi(k+ilk)-r(k+i)II~(;)+

L

lI.ó.û(k+ilk)II~(;) (2.11) ;=0

16CHAPTER 2. A STANDARD FORMULATION OF PREDICTNE CONTROL

Note that we do not necessarily start penalising deviations of z from r immediately, because there may be some delay between applying an input and seeing any effect. Note also that we assume that we know future rvalues exactly. This is reasonable in manyapplications. We will always assume that Hu S Hp, and that ~û(k

+

ilk) = 0

for i

2:

Hu, so that û(k

+

ilk) = û(k

+

Hu - 11k) for all i

2:

Hu. (Recall Figure 1.1.)

Most commonly we will assume that the controlled variables are the same as the measured outputs, so that Cz = Cy and Dz = Dy, and z(klk) = Y(klk).

In the following we use vec(O) to denote a vector, each element of which is O. We also use vec(a, b, c) to denote the vector [a, b, cjT, etc.

Constraints: E vec(~û(klk), ... ,~û(k

+

Hu - 11k), 1) S vec(O) F vec(û(klk), ... ,û(k

+

Hu - 11k), 1) S vec(O) G vec(.i(k

+

Hwlk), ... ,z(k

+

Hplk), 1) S vec(O) (2.12) (2.13) (2.14) Here E, F and G are matrices of suitable dimensions. We can use constraints of this form to represent possible actuator slew rates (2.12), actuator ranges (2.13),

and constraints on the controlled variables (2.14).

Example 2.1 Suppose we have a plant with 2 input variables and 2 controlled vari-ables. The control horizon is Hu

=

1 and the prediction horizon is Hp

=

2 (with Hw

=

IJ. We have the following constraints, which are to hold at each time:

-2 S ~Ul S 2 (2.15)

OS U2 S 3 (2.16)

Zl

2:

0 (2.17)

Z2

2:

0 (2.18)

3z1

+

5Z2 S 15 (2.19) To put these into our standard form, we rewrite the ~Ul constraints as follows.

-2 S

au, .. -au, s

2 ..

-~au,

-

1 S 0 .. [-1/2,0, -I[ [

~;:

1

s

0

Ju,

s

2 ..

~au,

-1

s

0 .. [I/2,0,-I[ [

~~:

1

s

0

- - - -- - -

-2.2. A STANDARD FORMULATION 17

so that we have

= [-1/2

0

-1]

E

1/2

0

-1

We leave finding the matrix F for this example to the reader (Problem 1), and go on to find G. We proceed in the same way, the main difference being that we have

to express the inequalities across the prediction horizon, which is greater than 1.

Z1 ~ 0 {} - Z1 :S 0

Z2 ~ 0 {} - z2 :S 0

1 1 [ "

1

3Z1

+

5Z2 :S 15 {} SZ1

+

"3Z2 - 1 :S 0 {} [1/5,1/3, -1]

7

from which we get

-1 0 0 0 0 [ ,,(H 11k) 0 0 -1 0 0 0 0 1/5 1/3 0 0 -1 z2(k

+

11k) 0 0 0 -1 0 0 z1(k

+

21k) :S 0 0 0 0 -1 0 z2(k

+

21k) 0 0 0 1/5 1/3 -1 1 0 v G :So

It will be seen that the dimensions of the matrices E, F and G can become very

large very quickly (Hp may be 10 or 20 or more). However, assembling them from

the given constraints is simpie, and this job can be automated. In fact, the software

which we will use - the Model Predictive Control Tooibox for Mat/ab-- does this

for the user.

The important thing about the constraints is that they are all in the form of linear

inequalities. When we come to solve the predictive control optimization problem, all

these inequalities will need to be translated into inequalities concerning ~û(k+ilk).

Since we assume a linear model, we will see that these inequalities remain linear,

even after this translation (Problem 2). This is going to be crucial for being able to

solve the optimization problem.

Notice that in the example there are variables (~U2,U1) which are not constrained.

It is also possible to have the converse situation, of variables which are constrained,

but do not appear in the cost function. This is only likely to occur with the Z

variables, and is represented in our standard formulation by having zero rows in

the weighting matrices Q(i). We therefore do not insist on Q(i)

>

0, but only on

18CHAPTER 2. A STANDARD FORMULATION OF PREDICTIVE CONTROL Feed Tank G w G P _Stock 1 -Headbox White water We Pap,

Figure 2.1: Paper machine with headbox

Example 2.2 Paper Machine Headbox

This example is based on an example in [30], and some of the software-based problems will refer to it.

Part of a paper-making machine is shown in Figure 2.1. The following varia bles are involved:

Inputs (u) States (x) Measured Controlled

Outputs (y) Variables (z) Stock ftowrateGp Feed tank level Hl

WW ftowrateGw Headbox level H2 H2 H2

Feed tank consistencyNI NI

headbox consistency N2 N2 N2

A linearised model of this machine has the (continuous-time) state-space matrices

o

o

-0.63 0.413

~

1

-

[

1.~74

o

Be -

1.34

-0.426 0

1.~741

-0.65

o

2.2. A STANDARD FORMULATION 19 where time has been measured in minutes. Predictive contral is to be applied to this plant, with an pdate interval of TB

=

2 minutes. Discretising the model with this sample interval (using the Matlab function c2d for example) gives

A

=

[~:~~!~ 0.4~66 ~ ~

1

B

=

[~:~:~~ ~:~:~~

1

o

0 0.2837 0 1.5237 -0.7391

0.1012 -0.6688 0.2893 0.4266 0.9929 0.1507

Note that the matrices Cv' Dy, Cz and Dz remain the same in the discrete-time model.

The prediction horizon is to be Hp

=

10 (that is, 20 minutes) and the control horizon

is to be Hu = 3. Since control oJ consistency is more important than control oJ level, errors in consistency (N2) are penalised more than errors in level (H2):

Q(i) =

[0

₀

2

n

Changes of the two valves are penalised equally, and the same pen alti es are applied at each step of the prediction and control horizons (that is, for each i):

R(i) = [0

0

6

0~6]

Since there is no delay in the plant (beyond the usuall-step delay - u(k) does not affect x(k)), errors are penalised /rom the first step, namely Hw = 1.

The valves are scaled so that their range oJ operation is:

(0 is usually the point at which the plant was linearised, so that the range is not always symmetrie about 0), and their slew rate (maximum rate oJ change) is such that they can move through 10% of the range in one step (2 minutes):

The headbox level (H2 ) is constrained, so that the headbox never runs dry, or

over-fills:

-3

:S

Zl

:S

5

Note that if an operating point is represented by u = 0, then with a linear model we must also have y = 0 at taht point (unless there is integration in the plant). Prom these constraints we can assem bIe the matrices E, F and G which appear in our standard Jormulation of the predictive control problem. These matrices will have the Jol/owing dimensions:

E : 4 x 7

(You should check this.)

20CHAPTER 2. A STANDARD FORMULATION OF PREDICTNE CONTROL

2.3 Prediction

In order to solve the predictive control problem, we must have a way of computing

the predicted values of the controlled variables, .î(k

+

ilk), from our best estimate of

the current state, x(klk), and the assumed fut ure inputs, or equivalently, the latest

input, u(k - 1), and the assumed future unput changes, t:.û(k

+

ilk).

It may seem that this is really part of the solution, not of the problem formulation, so

that it belongs in the next chapter. But it turns out that the way the predictions are

made has a great effect on the performance of the closed-loop system running under

predictive control. 80 the choice of prediction strategy is a 'tuning parameter' for

predictive control, just as choices of horizons and cost functions are. Furthermore,

the prediction strategy follows in a rather systematic way from assumptions made

about disturbanees acting on the system and measurement errors sueh as noise. 80

we can say that, rat her than ehoosing aprediction strategy, we are speeifying the

'environment' in whieh the plant is operating. And that properly belongs here, as part of the problem formulation.

However, prediction can get very complicated. 80 to avoid being too distracted by

it at this stage, we will deal with a few simple cases here - which will already cover

much of industrial practice - and come back to more general cases later in the

course.

2.3.1 No disturbances,

Full

state measurement

We will start with the simplest situation. Assume that the whole state vector is

measured, so that x(klk) = x(k) = y(k) (so Cy = I, Dy = 0). Also assume that we

know nothing ab out any disturbances or measurement noise. Then all we ean do is

to predict by iterating the model (2.8)-(2.10). 80 we get

x(k

+

11k)

=

Ax(k)

+

Bû(klk) x(k

+

21k)

=

Ax(k

+

11k)

+

Bû(k

+

11k)

=

A2x(k)

+

ABû(klk)

+

Bû(k

+

11k) x(k

+

Hplk)

=

Ax(k

+

Hp - 11k)

+

Bû(k

+

Hp - 11k) = AHpx(k)

+

AHp-l Bû(klk)

+ ... +

Bû(k

+

Hp - 11k) (2.20) (2.21 ) (2.22) (2.23) (2.24)

In the first line we have used û(klk) rather than u(k), because at the time when we

need to compute the predictions we do not yet know what u(k) will beo

Now recall that we have assumed that the input will only change at times k, k

+

1, ... , k

+

Hu - 1, and will remain constant af ter that. 80 we have û(k

+

ilk) =

~--~~---_ ...••...

_

... -_.- _. - - -

-2.3. PREDICTION 21

expressed in terms of !lû(k

+

ilk) rather than û(k

+

ilk), so let us do that now. Recall that !lû(k

+

ilk) = û(k

+

ilk) - û(k

+

i - 11k), and that at time k we already know u(k - 1). 80 we have

û(klk) = !lû(klk)

+

u(k - 1) û(k

+

11k) = !lû(k

+

11k)

+

!lû(klk)

+

u(k - 1) û(k

+

Hu - 11k) = !lû(k

+

Hu - 11k)

+ ... +

!lû(klk)

+

u(k - 1) and hence we get x(k

+

11k) = Ax(k)

+

B[!lû(klk)

+

u(k - 1)] x(k

+

21k) = A2x(k)

+

AB[!lû(klk)

+

u(k - 1)] B J!lû(k

+

11k)

+

!l~(klk)

+

u(k - 1)1, û(k+1lk) = A2x(k)

+

(A

+

I)B!lû(klk)

+

B!lû(k

+

11k)

+

(A

+

I)Bu(k - 1)

x(k

+

Hulk) = AHux(k)

+

(AHu-1

₊

_...

₊

_A

₊

_I)B!lû(klk)

₊

..

.

+

B!lû(k

+

Hu - 11k)

+

(AHu-l

+

..

.

+

A

+

I)Bu(k - 1)

Notice the change at this point

x(k

+

Hu

+

11k) = AHu+1x(k)

+

(AHu

+ ... +

A

+

I)B!lû(klk)

+

.

.. +

(A

+

I)B!lû(k

+

Hu - 11k)

+

(AHu

+

..

. +

A

+

I)Bu(k - 1)

x(k

+

Hplk) = AHpx(k)

+

(AHp-l

+ .

..

+

A

+

I)B!lû(klk)

+

22CHAPTER 2. A STANDARD FORMULATION OF PREDICTIVE CONTROL Finally we can write this in matrix-vector form:

x(k

+

11k) x(k

+

Hulk) x(k

+

Hu

+

11k) B AB+B A :AHu AHu+1 x(k)

+

o

o

B AB+B

..

future B

..

past [ t..û(klk)

1

(2.25) t..û(k

+

~u

- 11k)

If we assume that Dz

=

0 (which is the most common case) then the predictions of z are now obtained simply as

or

i(k

+

11k) = Czx(k

+

11k)

i(k

+

21k) = Czx(k

+

21k)

A warning is in order about the prediction equation (2.25). It involves computing Ai, possibly for quite large values of i. This can lead to numerical problems. If the plant is unstable, then some elements in Ai may become extremely large relative to others, and relative to elements in lower powers of A. Since computers work with finite-precision arithmetic, this can sometimes lead to wrong results. Similar problems can occur if the plant is stable; in this case some elements of Ai may become extremely small relative to others. Again wrong answers may result. (Using 'IEEE Standard' arithmetic, computers cannot distinguish between 1 and 1

+

I: if 11:1

<

10-16_{, approximately.)} (2.26) (2.27) (2.28) (2.29) (2.30)

2.3. PREDICTION 23

dek)

u(k)

~I

+

L _ _

PL_A_NT

_ _

:~~~O~

z(k)

1\ 1\

y(klk-l)

=

z(klk-l)

Figure 2.2: Output disturbance

The safest way of computing the predictions is probably to iterate one step at a time.

But this makes it difficult to fit the predictions into the optimization. Another

pos-sibility, which certainly reduces the number of flops (= 'fioating-point operations')

is to use 'Horner's method':

m

L

Ai

=

:t(A ... (A(A

+

I)

+

I) ... + I) +~ (2.32) i=O m times

=A(~Ai)

+I (2.33)

80 while the form of (2.25) is very convenient, and will be important in the next

chapter, it is not necessarily the form which should be used for implementation.

2.3.2 Constant output disturbanee

Now we will assume that there are disturbances acting on the plant. The simplest

assumption is that the measured and controlled outputs are the same (z = y), and

that there is an 'output disturbance' - see Figure 2.2:

y(k) = y(klk - 1)

+

d(klk)

=

Cyx(klk - 1)

+

d(klk)

(2.34) (2.35)

Another simple assumption is that this disturbance will continue unchanged durign

the prediction horizon. This means that, although we keep (2.9) and (2.10), we

predict y and z using

i(k

+

ilk)

=

y(k

+

ilk)

=

Cyx(k

+

ilk)

+

dek

+

ilk) (2.36)

where

dek

+

ilk) = d(klk) (2.37)

24CHAPTER 2. A STANDARD FORMULATION OF PREDICTIVE CONTROL

1. Measure the actual plant output y(k).

2. Estimate the disturbance as the difIerence between the actual and the estim-ated output.

3. Use that estimate to predict outputs over the prediction horizon.

We will always follow this scheme, but Step 2 will become more elaborate when we make more elaborate assumptions ab out the nature of the disturbance. In Step 3 we need to assume future input movements, in order to use (2.8) or (2.25) to predict x(k + ilk).

The assumption of a constant output disturbance, and the simplest disturbance estimation scheme of (2.36)-(2.37) is sometimes called the 'DMC scheme', because that is what is used in one of the leading prprietary predictive control products, DMC (= Dynamic Matrix Contraf) [12].

Notice that, even if Cy = I, we now have fi(k

+

ilk)

1'=

x(k

+

ilk). In general, this means that we need an observer to estimate 5;.

2.3.3 Using an observer

If we cannot measure the full state vector, or if the measured outputs consist of some linear combinations of the states, so that the states cannot be measured directly, then an observer can be used to estimate the state vector. It is instructive to see how the 'constant output disturbance' assumption we made in the previous subsection can be handled using observer theory. We can do this by 'augmenting' the model of the plant, so that it also includes a model of the output disturbance. We do this by defining a new 'state':

e(k) = [ x(k) ] d(k)

Since we assume that the output disturbance is constant, the new state and output equations are

[ d(k x(k

+

+

1)

1)]

= [A 0

0] [

I X(k)] d(k)

+

[ B ] 0 u(k) (2.38) [ x(k) ]

y(k) = [Cv

IJ

d(k) (2.39)

2.3. PREDICTION 25

---1

:

PLANT

1 1 y(k) u(k) 1 1 1 _ _ _ _

-OB SERVER

Figure 2.3: A State Observer

then applying the standard observer equation (see the Minitutorial) gives the fol-lowing estimates:

[ x(k d(k

+

+

1

11k) 1

k)] [(

= A 0

0)

I -

[Ll]

L2 [Cy IJ

] [

x(klk d(klk _

-1)] [

1)

+

B ]

0 u(k)

+

[ Ll ]

_L2 y(k) (2.41)

But in the previous subseetion we had

d(klk) = -Cyx(klk - 1)

+

y(k)

and sinee we assume that d(k

+

11k)

=

d(klk) we have

d(k

+ 11k)

= -Cyx(klk - 1)

+

y(k) (2.42) We ean see that the estimate obtained from the observer is the same as this if the observer gain matrix L is

This gives

[( A

_o

0) _ [

I

Ll ]

L2 [e Y

Il]

= [ -CA y

0]

0

The bloek-triangular strueture of this matrix shows that the observer's eigenvalues are those of the plant (the matrix A) and the remaining ones are zeros. The zeros

26CHAPTER 2. A STANDARD FORMULATION OF PREDICTNE CONTROL

Mini-Tutorial: State Observers

The general structure of a State Observer is shown in Figure 2.3, for a plant described by the equations

x(k

+

1)

=

Ax(k)

+

Bu(k), y(k) = Cx(k) (2.40)

It is a copy of the plant, with feedback from the measured plant output, through the gain matrix L, to correct the state estimate x.

The equations of the observer are:

x(k

+

11k) = Ax(klk - 1)

+

Bu(k)

+

L[y(k) - y(klk - 1)] y(klk - 1) = Cx(klk - 1)

Substituting the second equation into the first we get

x(k

+

11k) = (A - LC)x(klk - 1)

+

Bu(k)

+

Ly(k)

which is a stabie system if the eigenvalues of A - LC lie inside the unit disk. Furthermore if we define the state estimation error as e(k)

=

x(k) - x(klk -1), then using (2.8) we have

e(k

+

1)

=

(A - LC)e(k)

which shows that the state estimation error converges to zero if the observer is stabie, at a rate determined by the eigenvalues of A - LC.

If the pair (A, C) is observable, then given an arbitrary set of locations in the complex plane, a gain matrix L exists which places the observer eigenvalues at these locations. The problem of finding L is the dual of the state-feedback pole-placement problem, and can be solved using the same algorithm. (See functions place, acker, or kalman in Mat/ab's

Contro/ System Too/box.)

If the state and output equations of the plant are assumed to be subjected to white noise disturbances with known covariance matrices, then L can be chosen such that the mean square state estimation error is the smallest possible. The observer is then known as a Ka/man Filter.

2.3. PREDICTION 27 indicate that 'deadbeat' estimation of the disturbance is obtained, namely that the

disturbance estimate is exact after a finite number of steps - and we can see from

(2.42) that it is in fact exact af ter only 1 step (if the real disturbancebehaves exactly according to our model).

The fact that the observer's eigenvalues include those of A shows that this simple

disturbance estimation scheme can be used only with stable plants - otherwise the

estimate x(k+ 11k) will get worse and worse as k increases. 80 we see that, although

the 'DMC scheme' for estimating disturbances is simple and intuitive, it has some limitations. We also see that it is easy to overcome these limitations: we can keep

the same disturbanee model, but use a different observer gain matrix L. Even with

unstable plants, it is possible to find L such that all the observer's eigenvalues lie

inside the unit disk (providing that the pair (A, Cy ) is observable).

Note that if an observer is used, then the form of the state prediction equation (2.25)

remains the same, but x(k) has to be replaced by x(klk). We willlook at the details