Optimal Input Signal Design and MPC of Nonlinear Dynamical Systems

(1)

Optimal Input Signal Design and MPC of

Nonlinear Dynamical Systems

MARIETTE ANNERGREN

Masters’ Degree Project

Stockholm, Sweden September 2010

(2)

(3)

Abstract

The main topics of this master’s project are control theory, system iden-tification and convex optimization. The objective is to develop, imple-ment and test methods for optimal input signal design and for control

of a nonlinear dynamical system using MPC. The thesis begins with a

theoretical part, in which some known results in these fields are sum-marized. In the applied part of the thesis, methods are developed and exemplified in MATLAB.

Optimal input signal design is performed on specific FIR, ARX and DC-motor systems which are all controlled by an MPC. The implementation works very well for the FIR and ARX systems. The estimates of the true parameters fulfill the pre-specified requirements when using the optimal input signal constructed. An unexpected behavior is obtained of the estimates for the DC-motor system. Some additional approximations which were made in the design of the optimal input signal are thought to be the cause. Although, the source of the odd behavior were never confirmed. To be able to have a user-friendly environment for optimal input signal design, further work is necessary to overcome numerical problems in the implementation.

A general method of implementing MPC of nonlinear dynamical sys-tems is constructed. It is problematic to use MPC to control a nonlin-ear system. The reason for this is that a nonlinnonlin-ear system in general corresponds to a non-convex optimization problem in the MPC algo-rithm. Our method is based on making the problem convex through a linearization of the nonlinear system dynamics. The method is tested on simulations of a reaction wheel pendulum and a two link robot arm. It works very well and the systems fulfill the control objectives. Each optimization problem takes about 0.3-1 second to solve when using cvx. This is in some situations too slow to be able to control a system in reality. Further work is recommended on implementing our method with another solver so that it can be tested on an actual system which requires new updates every milli- or microsecond.

(4)

Design av optimal insignal och MPC för ickelinjära

dynamiska system

Detta examensarbete berör områden som reglerteknik, systemidentifier-ing och konvex optimersystemidentifier-ing. Syftet är att utveckla metoder för design av en optimal insignal och för att använda MPC på ickelinjära dynamiska system. Rapporten börjar med en teoretisk del som sammanfattar kän-da resultat inom dessa områden. Därefter följer beskrivningar av de metoder som har utvecklats och de exemplifieras genom simulering i MATLAB.

En metod för att designa en optimal insignal för specifika FIR-, ARX-och DC-system har utvecklats. Systemen styrs med hjälp av en MPC regulator. Metoden fungerar utmärkt för FIR- och ARX-system. Ett oväntat resultat erhålls för fallet med DC-motorn. Vi tror att det beror på de approximationer som har gjorts särskilt för detta system men det har inte kunnat bekräftas. För att skapa en användarvänlig miljö för design av en optimal insignal givet ett system och en regulator så krävs ytterligare arbete. Det bör fokusera på att eliminera de numeriska prob-lem som uppstår när metoden impprob-lementeras i MATLAB.

En allmän metod för att använda MPC på ickelinjära dynamiska sys-tem har implementerats. Ett ickelinjärt syssys-tem ger generellt upphov till ett ickekonvext optimeringsproblem i MPC algoritmen. Det är därför problematiskt att använda MPC för att reglera ickelinjära system. Vår metod bygger på att göra problemet konvext genom att linjärisera det ickelinjära systemet. Metoden har testats på simuleringar av en pendel och en två-länkad robotarm. Den presterar väldigt bra och systemen uppnår önskat beteende. Optimeringsproblemet tar cirka 0,3-1 sekund att lösa när vi använder cvx. Detta är i vissa fall för långsamt för att kunna reglera ett verkligt system. Ytterligare arbete där metoden imple-menteras med en annan lösare rekommenderas. Detta skulle möjliggöra att man kan testa den på ett verkligt system som kräver uppdateringar varje milli- eller mikrosekund.

(5)

3.1.5 The Relation Between System Identification- and Performance Degradation Cost . . . 15 3.2 Controllers . . . 16 3.2.1 PID-Controller . . . 16 3.2.2 MPC. . . 16 3.3 System Structures . . . 17 3.3.1 FIR System . . . 17 3.3.2 ARX System . . . 18 3.4 Mechanical Systems . . . 18 3.4.1 DC-motor . . . 18

3.4.2 Reaction Wheel Pendulum . . . 19

3.4.3 Two Link Robot Arm . . . 21

4 Analysis 23 4.1 Optimal Input Signal Design: V00 app(θ0) . . . 23

(6)

4.1.3 Minimum Volume Ellipsoid . . . 24

4.1.4 Accuracy . . . 24

4.2 Performance Degradation Cost: P-Control of FIR system . . . 25

4.2.1 FIR system . . . 26

4.2.2 Performance Degradation Cost . . . 27

4.3 Optimal Input Signal Design: MPC of FIR system . . . 27

4.3.1 FIR system . . . 27

4.3.2 MPC State Space Model . . . 28

4.3.3 Application Optimization Problem . . . 28

4.3.5 System Identification Cost. . . 29

4.3.6 System Identification Optimization Problem. . . 30

4.4 Optimal Input Signal Design: MPC of ARX system . . . 30

4.4.1 ARX System . . . 30

4.4.2 MPC State Space Model . . . 31

4.4.5 System Identification Cost. . . 33

4.4.7 Frequency Domain . . . 34

4.5 Optimal Input Design: MPC of DC-Motor System . . . 35

4.5.1 Discretized DC System Description . . . 35

4.5.2 Transfer Function . . . 36

4.5.5 System Identification Experiment. . . 37

4.6.1 Procedure . . . 38

4.6.2 Non-Convex Objective Function . . . 39

4.7 MPC of Reaction Wheel Pendulum . . . 40

4.7.1 Discretized Nonlinear System Description . . . 40

4.7.2 Linearized System Description . . . 40

4.7.3 Discretized Linearized System Description . . . 41

4.7.4 Objective Function . . . 41

4.8 MPC of Two Link Robot Arm. . . 42

4.8.1 Discretized Nonlinear System Description . . . 42

4.8.2 Discretized Linearized System Description . . . 42

4.8.3 Objective Function . . . 43

5 Examples 45 5.1 Optimal Input Signal Design: P Control of FIR system . . . 45

(7)

5.1.1 Example 1. . . 45

5.1.2 Example 2. . . 49

5.1.3 Example 3. . . 52

5.2 Optimal Input Signal Design: MPC of FIR system . . . 54

5.2.1 Example 1. . . 54

5.2.2 Example 2. . . 58

5.3 Optimal Input Signal Design: MPC of ARX system . . . 61

5.3.1 Example 1. . . 61

5.4 Optimal Input Signal Design: MPC of DC-motor . . . 64

5.4.1 Example 1. . . 64

5.4.2 Example 2. . . 67

5.5 MPC of Reaction Wheel Pendulum . . . 70

5.5.1 Example 1. . . 70

5.5.2 Example 2. . . 73

5.6 MPC of Two Link Robot Arm. . . 77

5.6.1 Example 1. . . 77

5.6.2 Example 2. . . 80

5.6.3 Example 3. . . 84

6 Discussion 89 6.1 Optimal Input Signal Design . . . 89

7 Conclusion 93 7.1 Optimal Input Signal Design . . . 93

(8)

(9)

Chapter 1

Introduction

This master’s project was carried out at the Department of Electrical Engineering, Stanford University, based on the research performed by Prof. Bo Wahlberg and Prof. Håkan Hjalmarsson at the Royal Institute of Technology and Ph.D. Candi-date for Prof. Stephen Boyd, Yang Wang, at Stanford University. The project is divided into two parts. The first part addresses the construction of optimal input signals for system identification given a model structure and a controller. This is called optimal input signal design. The second part considers control of nonlinear dynamical systems by implementing a so-called model predictive controller (MPC controller).

1.1 Main Topics

The main topics of this master’s project are control theory, system identification and convex optimization.

• “Control theory deals with analysis of dynamical systems and methodologies to construct controllers”, page 1 in [6]. A specific behavior of a dynamical system can be obtained by using a particular controller to regulate that system. One example is cruise control in cars where a controller ensures that the vehicle moves with constant speed. The main issues when controlling a system are to

– derive a mathematical model of the system, – design a controller.

In some situations it can be expensive and combined with risks to test con-trollers directly on the system. In those cases, computer simulations are used instead which requires a mathematical model of the system. Procedures in-volving system identification are used to construct those models.

(10)

• “System identification deals with the problem of building mathematical mod-els of dynamical systems based on observed data from the system”, page 1 in [10]. Dynamical systems can often be approximated by a mathematical model and by specific values of the parameters within that model structure. The main steps in the system identification process are to

– collect data from the system, – choose a structure for the model,

– choose a criterion which the model must fulfill to be accepted as a good approximation of the system,

– estimate values of the parameters in the model based on data, structure and criterion,

– validate the model.

The model is a good enough approximation of the system if it is validated. If not, a new model has to be constructed.

• Convex optimization is a useful tool for many practical applications such as control and system identification. Convex optimization problems belong to a certain class of mathematical optimization problems. They have a specific structure which make them qualify as convex. If a solution to a convex op-timization problem is found it is guaranteed to be globally optimal. Because of this, convex optimization problems can be solved reliably and efficiently, unlike optimization problems in general. One should always rewrite an opti-mization problem to convex form if possible, due to these qualities.

Prof. Stephen Boyd and Dr. Michael Grant at the Department of Electrical Engineering, Stanford University, have developed a programming framework for convex optimization problems called cvx, see [7]. cvx greatly simplifies the problem formulation compared to already existing optimization solvers. It is implemented in the programming tool MATLAB and uses MATLAB-syntax. cvxis used to solve the optimization problems in this master’s project.

1.2 Previous Work

1.2.1 Optimal Input Signal Design

The research performed by Prof. Bo Wahlberg and Prof. Håkan Hjalmarsson on optimal input signal design combines the objective of system identification with a pre-specified requirement on the system’s behavior. The estimated values of the parameters in the model should be such that the model can be validated in the system identification context and that the model fulfills the requirement set on the

(11)

1.3. OBJECTIVE

behavior in the control context. Roughly speaking, the goal of optimal input signal design should be to construct a model of a system which fulfills both the validation criterion and the requirement imposed on the control performance of the system. The benefit of taking the requested behavior into account already when modeling the system is that the number of possible models may be reduced, due to more constraints imposed on the model.

The underlying theory of optimal input signal design is well-developed but it exists no established framework for its implementation on specific modeling problems. It is desired to have a user-friendly programming environment which enables optimal input signal design and some examples of its implementation.

1.2.2 MPC of a Nonlinear Dynamical System

Different methods of applying MPC of a nonlinear dynamical system have been developed, see [11] and [5]. When implementing an MPC controller on a nonlin-ear system, one usually needs to solve a non-convex optimization problem in the MPC algorithm. In [11], this optimization problem is transformed into a convex one through a nonlinear change of variables. The obtained solution is thus guar-anteed to be globally optimal due to general properties of convex problems. It is extremely rare that such a transformation is possible. One other method is to solve the non-convex optimization problem as it is. Since we only know that the obtained solution is locally optimal, this will be an approximative solution. It is also quite difficult to construct a method of solving a non-convex optimization problem. The main issue for all methods is how to make the optimization problem solvable. There is no general algorithm which will transform non-convex optimization prob-lems into convex ones. It is desired to have a method of using MPC on a nonlinear dynamical system and it should be easy to implement.

1.3 Objective

The first objective of this master’s project is to construct a MATLAB-environment which designs optimal input signals for system identification given a model struc-ture and a controller. The strucstruc-tures are so-called FIR and ARX, and a specific mechanical model of a DC-motor. The MPC controller is used to regulate the systems. The second objective is to develop and document a general method of implementing MPC of nonlinear dynamical systems. The method is tested on two different simulated systems, a pendulum and a two link robot arm.

(12)

1.4 Acknowledgment

I would like to thank my supervisor Prof. Bo Wahlberg for his great support, invalu-able guidance and big inspiration throughout this master’s project. I also wish to thank my second supervisor Ph.D. Candidate Yang Wang, who with a great deal of patience helped me with good advice, numerous guiding discussions and an unyield-ing positive attitude. The main part of this master’s project was carried out at the Stanford University, for this rare and much appreciated opportunity I would like to thank Prof. Stephen Boyd. My time at Stanford University has been a tremendous learning experience for me. Lastly, I would like to thank Denise Murphy for her warm welcome and for making sure that I had a place to work and a computer to use at the Department of Electrical Engineering throughout my time there.

(13)

Chapter 2

Background

2.1 Optimal Input Signal Design

The objective is to construct an optimal input signal. The input signal is optimal in the sense that, when used in a system identification experiment to estimate the values of the parameters in the system, it ensures with a given probability that the resulting model fulfills the requirement on the system’s behavior in the control con-text. The fulfillment of the requirement should ensure that the model has the same behavior as the system. It is not necessary to achieve all the characteristics of the system. A choice should be made regarding what qualitative behavior is important and what may be neglected. The requirement depends on the purpose of the system and thus may vary from application to application.

The focus of this master’s project is to construct optimal input signals for system identification. The questions addressed are:

• How can the characteristics of an optimal input signal be determined? • How can the estimated model be constructed?

A theoretical background of optimal input signal design is presented in Chapter3. The methods developed and implemented to perform optimal input signal design are described in Chapter4. A set of system identification experiments was performed during this master’s project to exemplify some of the aspects of optimal input signal design. The descriptions and results of these are given in Chapter5.

2.2 MPC of Nonlinear Dynamical Systems

One of the benefits with MPC is that it does not require much physical knowledge about the system to be able to regulate it. It is sufficient to only have a math-ematical model of the system. Other controllers, for instance the PID-controller, demands a great deal of knowledge of the system’s behavior. In PID control design,

(14)

the user adjusts the parameters to handle control objectives and system constraints indirectly. In MPC, the user specifies these explicitly as part of the optimization problem. Thus, MPC requires less intuition and experience of its user to get good control performance. The fact that MPC only requires a mathematical model en-ables a generalization to be made of the implementation method, meaning that the same method can be used on different systems. It is therefore desirable to find a procedure for implementing MPC of a nonlinear dynamical system since this would lead to a general method for several nonlinear systems.

The questions addressed are:

• Which approximations need to be made in the model of the nonlinear dynam-ical system to enable use of an MPC controller?

• How do these approximations affect the behavior and level of control of the nonlinear dynamical system?

A theoretical background of MPC is presented in Chapter 3. The method of im-plementation used in this master’s project is described in detail in Chapter4. The method is exemplified by two control experiments. Descriptions and results of these experiments are given in Chapter5.

(15)

Chapter 3

Theory

This chapter describes existing theory which this master’s project is based on. It is necessary to have some knowledge about system identification and control theory to be able to follow the text in the following sections.

3.1 Optimal Input Signal Design

The theory described in this section is a summary of both standard results on system identification and more recent work on optimal input signal design. For a full description of the underlying theory see [10] and [8].

3.1.1 Notation

The notation used follows [10] and [8]. Here is a short list of the symbols: AsN() = asymptotic normal distribution

d = dimensions of the parameter vector θ

DM = set of values over which θ ranges in a model structure

Ex = mathematical expectation of a random vector x

¯

Ex(t) = limN →∞_N1 PNt=1Ex(t)

M = model structure (a mapping from a parameter space to a set of models)

N = number of measurements performed on the system Pθ = asymptotic covariance matrix of θ

q, q−1 = forward and backward shift operators u(t) = input variable at time t (input signal) VN(θ, ZN) = criterion function to be minimized

¯

V(θ) = limit of criterion function as N → ∞ v(t) = disturbance variable at time t

y(t) = output variable at time t (output signal)

ˆy(t|θ) = predicted output at time t using model M(θ) and based on Zt−1

(16)

ZN = {y(0), u(0), y(1), u(1) . . . , y(N), u(N)} (t, θ) = y(t) − ˆy(t|θ), prediction error

λ = used to denote variance

θ = vector used to parametrize models θ0 = true parameter values of a given system

ˆθ_N _{= estimate of θ based on Z}N

θ∗ = limiting estimate of θ ϕ(t) = regression vector at time t

All vectors are assumed to be column vectors if they are not defined otherwise. The gradient of a function f(x) is denoted f0_{(x) and is a column vector. The hessian of}

the same function is denoted f00_(x).

3.1.2 Requirements on a Model

The model of a system should have a behavior that is sufficiently close to the be-havior of the system. There are different methods for determining what sufficiently close means. One of them is based on minimizing the prediction error. This method is used throughout this master’s project and is described in detail in the following section.

Prediction Error

A set of data, denoted ZN_{, is collected from the system by performing measurements}

and recording the results. This is all done in discrete-time. The data set consists of the measured output signal y(t) from the system corresponding to a particular input signal u(t) at different times, denoted t. If the same sequence of input signals is fed into a candidate model instead of the system itself, a set of predicted output signals ˆy(t|θ) is obtained. In other words, the output signals from the candidate model predict the output signals from the system. The symbol θ denotes the vector of parameters in the system that are to be estimated. Given a particular candidate model with θ = θ∗, the prediction error is defined as

(t, θ∗) = y(t) − ˆy(t|θ∗). (3.1)

A good candidate model should give a small prediction error. On page 198 in [10], the principle of parameter estimation is phrased as:

“Based on Zt_{we can compute the prediction error (t, θ) using [Equation (3.1)]. At}

time t = N, select ˆθN so that the prediction errors (t, ˆθN), t = 1, 2, ..., N, become

as small as possible.”

This means that the sequence of prediction errors should be minimized with respect to θ to be able to obtain a good candidate model. A qualitative measurement of the size of the prediction errors is necessary to be able to perform the minimization.

(17)

3.1. OPTIMAL INPUT SIGNAL DESIGN

Thus, it is desired to have a mapping from the prediction error sequence to a scalar value. Such a function is called a criterion function and is denoted VN(θ, ZN). This

leads to the following definition of ˆθN:

ˆθ_N _{= ˆθ}_N_(ZN_{) = arg min V}

N(θ, ZN), θ ∈ DM. (3.2)

There exist several functions that can be used as the criterion function VN(θ, ZN).

The function used here is the least-squares criterion. Least-Squares Criterion

Let the candidate models have a linear regression structure, meaning that the func-tion that calculates the predicted output signals is a scalar product between a known data vector ϕ(t) and the parameter vector θ, see page 8 in [10],

ˆy(t|θ) = ϕT_(t)θ. _(3.3)

The prediction error, defined in Equation (3.1), then becomes

(t, θ) = y(t) − ϕT(t)θ. (3.4)

The least-squares criterion for the linear regression structure in Equation (3.3) is defined as VN(θ, ZN) = 1 N N X t=1 (y(t) − ϕT_(t)θ)2_. _(3.5)

Convergence Properties of the Least-Squares Estimate

The minimization problem in Equation (3.2) can be solved analytically when the criterion function is defined as in Equation (3.5). By differentiating VN(θ, ZN) with

respect to θ and setting the obtained expression equal to zero, one gets ˆθLS N = arg min VN(θ, ZN) = " 1 N N X t=1 ϕ(t)ϕT(t) #−1 1 N N X t=1 ϕ(t)y(t), (3.6) with ˆθLS

N denoting the least-squares estimate. Now let

R(N) = 1 N N X t=1 ϕ(t)ϕT(t) (3.7)

(18)

and f(N) = 1 N N X t=1 ϕ(t)y(t), (3.8)

where R(N) is a d-by-d matrix and f(N) is a column vector of length d. Assume that the output signals from the system are generated by the linear regression structure

y(t) = ϕT(t)θ0 + v0(t). We use the convention that a subscript “0” denotes the

true value of a parameter, so θ0 is the true value of θ. Similarly, v0 denotes the

disturbance signal that actually is affecting the system. It is added to the structure since no real system is undistorted by noise. Thus, the estimate in Equation (3.6) becomes ˆθLS N = [R(N)] −1 1 N N X t=1 ϕ(t) ϕT(t)θ0+ v0(t) = θ0+ [R(N)]−1 1 N N X t=1 ϕ(t)v0(t). (3.9)

The estimate ˆθN should be close to the true value θ0 and it should converge to this

value as more data samples are collected, meaning as N → ∞. If the sequence v0(t)

is small in comparison to ϕ(t) the estimate will be close to θ0. This statement is

easily proved by substituting Equation (3.7) into Equation (3.9). Certain constraints must be put on the error term

[R(N)]−1 1 N N X t=1 ϕ(t)v0(t)

for the estimate to converge to the true value. Sufficient conditions are stated on page 205 in [10] and are

• R(N) → R∗ _{= ¯}_Eϕ_(t)ϕT_{(t) as N → ∞, R}∗ _{should be nonsingular,}

• 1

N

PN

t−1ϕ(t)v0(t) → f∗= ¯Eϕ(t)v0(t) as N → ∞, f∗ should be zero.

Thus, if these constraints are met ˆθLS

N will converge to θ0, that is

ˆθLS

(19)

Statistical Properties of the Least-Squares Estimate By the definition in Equation (3.2), all estimates ˆθN should fulfill

V_N0(ˆθN, ZN) = 0. (3.11)

Assume that θ∗ _{is a unique global minimizer of ¯}_V_{(θ), which is the limit of the}

criterion function VN(θ, ZN) as N → ∞. A Taylor expansion of Equation (3.11)

around the point θ∗ _yields

0 = V0

N(θ

∗_{, Z}N_{) + V}00

N(ξN, ZN)(ˆθN − θ∗), (3.12)

where ξN is a vector whose values lie between those of ˆθN and θ∗. It can be shown

that V00

N(ξN, ZN) → ¯V00(θ∗) as N → ∞ with probability 1, see page 281 in [10].

Thus, the Taylor expansion for large N can be rewritten as (ˆθN − θ∗) = −[VN00(θ ∗_{, Z}N_)]−1_V0 N(θ ∗_{, Z}N_). _(3.13) The vector −V0 N(θ∗, ZN) is given by − V_N0 (θ∗, ZN) = 2 N N X t=1 ψ(t, θ∗)(t, θ∗), (3.14) where ψ(t, θ∗) = −d dθ(t, θ)|θ=θ∗. (3.15)

Since θ∗ _{is a unique global minimizer of ¯}_V_{(θ), one has that}

¯

V0(θ∗) = 2 ¯Eψ(t, θ∗)(t, θ∗) = 0. (3.16)

Assume that Equation (3.14) is a sum of random variables ψ(t, θ∗_{)(t, θ}∗_{) with zero}

mean value and that V0

N(θ∗, ZN) has an increasing independence between distant

terms with respect to time. Then under certain conditions, which we assume to be true, the well-known central limit theorem implies that

r 2 N N X t=1 ψ(t, θ∗)(t, θ∗) ∈ AsN(0, Q), (3.17) where

(20)

Q= lim

N →∞

N

2 · E{[VN0(θ∗, ZN)][VN0 (θ∗, ZN)]T}, (3.18)

see page 282 in [10]. This means that q

2

N

PN

t=1ψ(t, θ∗)(t, θ∗) converges in

distri-bution to the normal distridistri-bution with zero mean and covariance matrix Q. Based on Equations (3.13), (3.14), (3.17) and (3.18), one obtains

s

N

2(ˆθN − θ∗) ∈ AsN(0, Pθ), (3.19)

where

Pθ = [ ¯V(θ∗)00]−1Q[ ¯V(θ∗)00]−1, (3.20)

see page 282 in [10]. This means thatqN

2(ˆθN− θ

∗_{) converges in distribution to the}

normal distribution with zero mean and covariance matrix Pθ.

The least-squares estimate ˆθLS

N has Pθ = λ0[2R∗]−1, see page 284 in [10], where

λ0 is the variance of v0. The expression for Pθ is obtained by a straightforward

calculation based on Equation (3.14) and Equation (3.18). Thus, s

N

2 (ˆθLSN − θ∗) ∈ AsN(0, λ0[2R∗]−1). (3.21)

The χ2_{-distribution can be used to describe the convergence instead. A normally}

distributed variable with zero mean and variance σ2 _{will become χ}2_{-distributed if}

squared and divided by σ2_{, with the degrees of freedom equal to the dimension of}

x. Formally, if x ∼ N(0, σ2) then x_σ22 ∼ χ2(d) where d = dim(x). For the situation

in Equation (3.21) this means

N

2λ0(ˆθ LS

N − θ∗)T2R∗(ˆθLSN − θ∗) ∈ Asχ2(d). (3.22)

3.1.3 Requirements on Application Performance

When creating a model of a true system it is necessary that the specifications on the application are met by the model. The estimates of the unknown parameters θ of the system must lie within a range of their actual values which ensures that the model has the same qualities as the true system. The function used to confirm if the model meets the specifications is called performance degradation cost. There are different ways to express this function depending on which qualities of the system one would like to check.

(21)

3.1.4 Geometric Approach

The results obtained from the system identification theory can be demonstrated in a geometric framework using confidence ellipsoids. The relation between the requirements on the model and the requirements on the application performance are stated and explained in the following sections.

System Identification Cost

The so-called system identification cost VSI(θ) should have the following qualities:

• V0

SI(θ0) = 0,

• V00

SI(θ0) 0.

The values of the estimated parameters should be chosen such that VSI(θ) is

mini-mized, meaning ˆθN = arg minVSI(θ). The same convergence of ˆθN as in Equation

(3.22) is obtained if VSI(θ) is chosen such that VSI00(θ) = 2R∗. These ˆθN will converge

in distribution as

N

2λ0(ˆθN

− θ₀)TV_SI00(θ0)(ˆθN − θ0) ∈ Asχ2α(d), (3.23)

where α corresponds to the α-percentile of the χ2_{-distribution and d to its degrees}

of freedom, i.e. the number of estimated parameters, see Equation (3.22). This convergence and the other specified qualities are achieved by choosing VSI(θ) =

¯

V(θ) − λ0, with VN(θ, ZN) defined as in Equation (3.5).

Second Order Approximation

A second order approximation can be made of the system identification cost by making a Taylor expansion of VSI(θ) around θ0. The approximation becomes

VSI(θ) ≈ VSI(θ0) + VSI0 T(θ0)(θ − θ0) +

1

2(θ − θ0)TVSI00(θ0)(θ − θ0) + O((θ − θ0)3). (3.24)

Since VSI(θ0) = 0 and V_SI0 (θ0) = 0, the Taylor expansion in Equation (3.24) can be

rewritten as

VSI(θ) ≈ 1₂(θ − θ0)TVSI00(θ0)(θ − θ0). (3.25)

One can see that the second order approximation in Equation (3.25) corresponds to a scaled left hand side of Statement (3.23). Hence, the approximation of the system identification cost in Equation (3.25) will belong to the χ2_{-distribution.}

(22)

Ellipsoid and Contour Plots

A geometric interpretation of the system identification cost can be made. The set of admissible θ for a certain threshold forms an ellipsoid,

E_V_SI = {θ | 1 2λ0(θ − θ0) T_V00 SI(θ0)(θ − θ0) ≤ χ2_α(d) N }, (3.26)

see page 558 in [10]. The ellipsoid EVSI is a confidence ellipsoid. The statistical

in-terpretation of this is that if the least-square estimate of the parameters, ˆθLS

N , were

to be made a number of times and the confidence ellipsoid constructed for each of those times, the obtained regions would contain the true values of the parameters,

θ0, with probability α. Another statistical interpretation is that the estimate of θ,

obtained when using the least-squares criterion, will lie inside of EVSI with

proba-bility α.

The directions of the semi-axis of the ellipsoid are determined by the eigenvectors of V00

SI(θ0), their lengths by the corresponding eigenvalues. Note that VSI00(θ0) is

a symmetric matrix, thus its eigenvalues will be real and its eigenvectors will be orthogonal to each other.

Performance Degradation Cost

The performance degradation cost, denoted Vapp(θ), should have the following

prop-erties: • Vapp(θ0) = 0. • V0 app(θ0) = 0, • V00 app(θ0) 0.

The maximum allowed value of the performance degradation cost is limited by the required accuracy of the model, defined as 1/γ.

Second Order Approximation

A second order approximation can be made of the performance degradation cost by use of Taylor expansion of Vapp(θ) around θ0. The approximation becomes

Vapp(θ) ≈ Vapp(θ0) + Vapp0 T(θ0)(θ − θ0) +

1

2(θ − θ0)TVapp00 (θ0)(θ − θ0) + O((θ − θ0)3). (3.27)

(23)

Vapp(θ) ≈ 1

2(θ − θ0)TVapp00 (θ0)(θ − θ0). (3.28)

Ellipsoid and Contour Plots

A geometric interpretation can also be made of the performance degradation cost. The set of admissible θ for a certain accuracy forms an ellipsoid defined as

EVapp = {θ |

1

2(θ − θ0)TVapp00 (θ0)(θ − θ0) ≤ 1

γ}. (3.29)

3.1.5 The Relation Between System Identification- and Performance

Degradation Cost

System identification methods are used to obtain as good estimates as possible of the true parameters of the system, but the requirements set on the true system’s behavior should also be met by the model. The second order approximation of the system identification cost in Equation (3.25) and the distribution in Equation (3.23) imply

VSI(θ) ≤

λ0χ2α(d)

N , (3.30)

and the second order approximation of the performance degradation cost in Equa-tion (3.28) along with the required accuracy imply

Vapp(θ) ≤ 1

γ. (3.31)

It is desirable to have a criterion that ensures that whenever Equation (3.30) holds, Equation (3.31) holds, too. In other words, when a good candidate model has been found it will also fulfill the requirements set on its behavior. The geometric interpretation of this statement is that

EVSI ⊆ EVapp. (3.32)

It is shown in [8] that the inequality in Equation (3.32) holds if and only if

VSI(θ) ≥

γλ0χ2α(d)

N Vapp(θ). (3.33)

Another condition which ensures that the inequality in Equation (3.32) holds can be derived from the expressions of the parametric regions EVSI and EVapp. The

(24)

V_SI00(θ) γλ0χ 2 α(d) N V 00 app(θ). (3.34)

3.2 Controllers

Controllers are a tool for regulating dynamical systems so that desirable behavior is obtained. The goal could be to create an output signal from the system that has a value which is close to a particular set point or to minimize the effect of disturbances on the system. There exist several different types of controllers. The PID- and MPC controller will be discussed in detail in the following sections.

3.2.1 PID-Controller

PID stands for proportional-integral-derivative. The PID-controller uses the differ-ence between the output signal and the set point, the so called error, to create an appropriate input signal to the system. The input signal should make the system compensate for the existing error and hence make it smaller. The PID-controller consists of three parts given by its name, one proportional, one integral and one

derivative. It can be described by

u(t) = KPe(t) + KI

Z t

t0

e(τ)dτ + KDde(t)_dt , (3.35)

where u(t) is the input signal to the system, e(t) is the error and KP, KI and KD

are the constants which scale the influence of each part of the controller respectively. The main issue is to tune these constants so that the best possible behavior of the system is obtained. The system should be stable and have an output signal close to the set point.

A controller consisting of only the proportional part is used in this master’s project, i.e. a P-controller.

3.2.2 MPC

MPC is based on

• a model of the system,

• measured data from the system,

• a cost function which penalizes undesirable behavior, • constraints which represent physical system limits.

(25)

3.3. SYSTEM STRUCTURES

One can predict future output signals from the system based on the current mea-sured data and the mathematical model. These predictions can be described as a function of a future input signal sequence implemented on the system. They are used in the cost function. The cost function is minimized with respect to the se-quence of input signals. The first input signal of that sese-quence is applied to the system. New measured data is obtained and the procedure of calculating the ap-propriate input signal starts over again. The steps of MPC can be summarized as

1. obtain measured data from the system,

2. use current data and the model of the system to predict future output signals as a function of future input signals,

3. insert the predictions in the cost function and minimize it with respect to a future input signal sequence,

4. apply the first input signal in the obtained optimal input signal sequence, 5. repeat procedure from step 1.

The prediction horizon is how far ahead the MPC algorithm will predict output signals. It is denoted T .

3.3 System Structures

The first step when performing system identification is to decide which kind of structure the model of the system should have. The choice is often based on the physical characteristics of the system, empirical results and the experience of the user. Two different kinds of systems will be described in this section.

3.3.1 FIR System

FIR stands for finite impulse response. The system can be described as

y(t) = b1u(t − 1) + b2u(t − 2) + ... + bnbu(t − nb) + e(t), (3.36)

where y(t) is the output signal, u(t) is the input signal and {e(t)} is white noise with variance λ.

The system is called finite impulse response because its response to an impulse will be zero after a finite number of time steps.

(26)

3.3.2 ARX System

The ARX system can be described as

y(t) − a1y(t − 1) − a2y(t − 2) − ... − anay(t − na) =

b1u(t − 1) + b2u(t − 2) + ... + bnbu(t − nb) + e(t). (3.37)

The FIR system defined in Equation (3.36) is a special case of Equation (3.37) with

na= 0.

3.4 Mechanical Systems

Simulations are made of mechanical systems to be able to exemplify the theory re-garding both optimal input signal design and implementation of MPC on a nonlinear dynamical system. The mechanical systems used are a DC-motor, a pendulum and a two link robot arm. The mathematical model of these are explained in detail in the following sections.

3.4.1 DC-motor

The model of the DC-motor used comes from [4].

s



_R

, J

_R



_1,

J

₁



_2,

J

₂

L

_a

R

_a

u

s

Figure 3.1. Diagram over DC-motor system.

The diagram of the DC-motor is displayed in Figure3.1. The notation in the figure means

(27)

3.4. MECHANICAL SYSTEMS

u = voltage,

Ra = resistance in armature,

La = inductance of armature,

i = current,

s = back electromotive force, φi = rotation angle of wheel i,

Ji = moment of inertia of wheel i,

for i = 1, 2, R.

State Space Model

The system description given on pages 15-16 in [4] is used here. The states of the system are x(t) = [φ1(t) ˙φ1(t) φ2(t) ˙φ2(t)]T∈ R4. The system dynamics is given by

˙x(t) =       0 1 0 0 −kS P1 −RaP2−f_T2kakv RaP1 kS P1 fS P1 0 0 0 1 kS J2 fS J2 −kS J2 −P3 J 2       x(t) +      0 fTka RaP1 0 0      u(t), y(t) = [0 0 1 0]x(t). (3.38) where

fi = viscous friction of wheel i,

fS = damping coefficient of spring,

fT = transformation factor from rotor to first wheel,

ka = rotor torque constant,

kS = elasticity of spring,

kv = back electromotive force constant,

P1 = J1+ JRfT2,

P2 = fRfT2 + f1+ fS,

P3 = fS+ f2,

for i = 1, 2, R.

3.4.2 Reaction Wheel Pendulum

A pendulum is connected to a hanging device in one end. On the other end there is a reaction wheel by which the movement of the pendulum can be controlled. The model used here is based on [1].

(28)

p

_r

m_r m_p

Figure 3.2. Reaction wheel pendulum.

The reaction wheel pendulum is displayed in Figure3.2. The notation in the figure means

mp = mass of pendulum,

mr = mass of reaction wheel,

θp = angle of pendulum,

θr = angle of reaction wheel.

Nonlinear System Description

The nonlinear system description given in [1] is used here. The states of the system are x(t) = [θp(t) ˙θp(t) θr(t) ˙θr(t)]T ∈ R4. We will also use the vectors xp(t) =

[θp(t) ˙θp(t)]T ∈ R2, xr(t) = [θr(t) ˙θr(t)]T ∈ R2 and u(t) ∈ R. The system dynamics

are described by ¨θp(t) = − mgl J sin(θp(t)) − k Ju(t) (3.39) and ¨θr(t) = k Jr u(t), (3.40) where

(29)

3.4. MECHANICAL SYSTEMS

J = Jp+ mpl2p+ mrlr2,

Jr = moment of inertia of the reaction wheel around its center of mass,

Jp = moment of inertia of the pendulum around its center of mass,

k = motor torque constant,

l = distance from pivot to the center of mass of pendulum and reaction wheel, lp = distance from pivot to the center of mass of the pendulum,

lr = distance from pivot to the center of mass of the reaction wheel,

m = mass of pendulum and reaction wheel.

3.4.3 Two Link Robot Arm

A robot arm is attached to a hanging device in one end. It has two motors located at each of its two links. Torque can be applied to the system by using the motors. The robot arm is moving in a horizontal plane, meaning that any impact from gravitation on the system can be disregarded. The model used is partly based on [9]. l_2,m₂ ₁ ₂ l_1,m₁ 1 2

Figure 3.3. Robot arm.

(30)

li = length of link i,

mi = mass of link i,

θi = angle of link i,

τi = input signal (torque) applied on link i,

for i = 1, 2.

Nonlinear System Description

The nonlinear system description is derived in [9]. The states of the system are

x(t) = [θ1(t) ˙θ1(t) θ2(t) ˙θ2(t)]T ∈ R4. We will also use the vectors τ(t) = [τ1(t) τ2(t)]T

∈ R2 _{and Θ = [θ}

1(t) θ2(t)]T∈ R2. The dynamics of the robot arm can be described

as τ1(t) = (m1+ m2)l21¨θ1(t) + m2l22(¨θ1(t) + ¨θ2(t)) + m2l1l2(2¨θ1(t) +¨θ2(t))cos(θ2(t)) − m2l1l2(2 ˙θ1(t) ˙θ2(t) + ˙θ22(t))sin(θ2(t)) (3.41) and τ2(t) = m2l22(¨θ1(t) + ¨θ2(t)) + m2l1l2¨θ1(t)cos(θ2(t)) −m2l1l2˙θ1(t)2sin(θ2(t)). (3.42)

The dynamic Equations (3.41)-(3.42) can be rewritten in a matrix form as

M(Θ(t)) ¨Θ(t) + W (Θ(t), ˙Θ(t)) ˙Θ(t) = τ(t), (3.43) with M(Θ(t)) = (m1+ m2)l 2 1+ m2l22+ 2m2l1l2c2(t) m2l22+ m2l1l2c2(t) m2l22+ m2l1l2c2(t) m2l22 ! (3.44) and W(Θ(t), ˙Θ(t)) = −2m2l1l2˙θ2(t)s2(t) −m2l1l2˙θ2(t)s2(t) −m2l1l2˙θ1(t)s2(t) 0 ! , (3.45)

(31)

Chapter 4

Analysis

Here follows a description of the different methods and approaches which are con-sidered and then implemented to be able to reach the objectives of this master’s project.

4.1 Optimal Input Signal Design: V

_app00

(θ

0

)

There exist different methods for determining the hessian of the performance degra-dation cost. Three such methods have been implemented in this master’s project. They are denoted numerical hessian, quadratic approximation and minimum

vol-ume ellipsoid. All three methods require descent estimates of the true parameters θ0 to be able to construct the hessian of the performance degradation cost. The

descent estimates are simply denoted θ0 in the following sections.

4.1.1 Numerical Hessian

The numerical hessian is used to approximate V00

app(θ0). Different values of Vapp(θ)

are obtained by simulating the model for different values of θ. These values of

Vapp(θ) are then used in the well-known central difference equation to approximate

the hessian. The limiting factor of this approach is that no general way was found of constructing the numerical hessian for each order of the system, i.e. for each number of unknown parameters. The expression for the hessian depends significantly on the number of parameters to be estimated in the system. This leads to a lot of time consuming work were each possible case demands a uniquely defined expression for the hessian.

4.1.2 Quadratic Approximation

The expression for V00

app(θ) is approximated by a function which is quadratic in θ.

(32)

fapprox(θk) = θTkQθk+ 2qTθk− θT0Qθ0−2qTθ0. (4.1)

The constant matrix Q and the constant vector q are determined by minimizing the difference between the value of Equation (4.1) and the value of the true function

Vapp(θ) for a set of points θk surrounding the true point θ0. This is equivalent to

the optimization problem minimize

K

X

k=1

(fapprox(θk) − Vapp(θk))2, (4.2)

with variables Q and q. The approximation of the hessian of Vapp(θ) then becomes

2Q, which is the second derivative of fapprox(θk) with respect to θk. One difficulty

with this approach is that it is not clear how to choose the set of points θk.

4.1.3 Minimum Volume Ellipsoid

A set of points θk which fulfills the criterion on the performance degradation cost,

see Equation (3.29), and which surrounds the true point θ0 is chosen. The minimal

ellipsoid containing all these points is constructed and is then used as an approxi-mation of the region EVapp. The calculation of the minimal ellipsoid containing all

points θk is equivalent to the optimization problem

maximize log det A w.r.t. A and b, s.t. ||AθT

k + b||22 ≤1 for all k, (4.3)

where (det A−1₎1/2 _{is proportional to the volume of the ellipsoid. We can therefore}

maximize the volume of the ellipsoid by maximizing det A. This is a convex opti-mization problem and thus solvable. The approximation of the hessian of Vapp(θ)

then becomes 2AT_A/γ_{. For more details about Optimization Problem (4.3) see}

pages 410-411 in [3]. One difficulty with this approach is, as in the quadratic ap-proximation, that it is not clear how to choose the set of points θk.

4.1.4 Accuracy

Two methods were used to check the accuracy and reliability of the three approaches described above. The first method plots the contours of Vapp(θ) for all values of θ

such that Vapp(θ) = 1/γ. A comparison could then be made of the true parameter

region for which one has an acceptable performance degradation cost, and the ap-proximated ellipsoidal region.

In the second method a set of points θk surrounding the true point θ0 is chosen.

(33)

4.2. PERFORMANCE DEGRADATION COST: P-CONTROL OF FIR SYSTEM

1/γ − < Vapp(θk) < 1/γ + , (4.4)

for some fixed > 0. Thus, every point θk is either classified as inside the

inter-val defined by the inequalities in Equation (4.4) or outside it. The three different approximations of the hessian are then tested by calculating the fraction of points that get classified in the same way by both the true Vapp(θ) and the second order

approximation.

After testing all three methods, we found that the numerical hessian gives the best approximation. This is why this method is used in the examples in Chapter 5.

4.2 Performance Degradation Cost: P-Control of FIR

system

In the examples in Section5.1the optimal input signal is designed for an FIR system controlled by a P controller.

Figure 4.1. P-controller implemented on FIR system.

The complete system which is simulated is displayed in Figure 4.1. The notation used in the figure means

(34)

d = disturbance signal, r = reference signal, u = input signal,

y = measured system output signal, z = system output signal.

4.2.1 FIR system

The model structure is described by the linear difference equation

y(t) = b1u(t − 1) + b2u(t − 2) + · · · + bnbu(t − nb) + e(t). (4.5)

The adjustable parameters are

θ= [b1 b2. . . bnb]

T _(4.6)

and the regression vector is

ϕ(t) = [u(t − 1) u(t − 2) . . . u(t − nb)]T. (4.7)

The transfer function from the input signal to output signal, i.e. from u(t) to z(t), is

G(q) = b1q−1+ b2q−2+ · · · + bnbq

−nb_, (4.8)

when no disturbance d(t) is added to the system. A P controller with a constant value K is used in a negative feedback loop, i.e. u(t) = −K(y(t) − r(t)), leading to the following expression for the closed loop output signal:

z(t) = −K(b1q−1+ b2q−2+ · · · + bnbq

−nb)(y(t) − r(t)). (4.9)

The closed loop transfer function, i.e. the transfer function from the reference signal to the output signal, is

G(q) = K(b1q −1_{+ b} 2q−2+ · · · + b−nbq −nb 1 + K(b1q−1+ b2q−2+ · · · + bnbq −nb. (4.10)

when no disturbance d(t) is added to the system.

A disturbance is added to the system according to Figure4.1. Thus, Equation (4.9) is modified to

y(t) = −K(b1q−1+ b2q−2+ · · · + bnbq

(35)

4.3. OPTIMAL INPUT SIGNAL DESIGN: MPC OF FIR SYSTEM

4.2.2 Performance Degradation Cost

The performance degradation cost Vapp(θ) is defined as

Vapp(θ) = lim t→∞ y(θ, t) − y(θ0, t) y(θ0, t) !2 . (4.12)

This can be called a relative performance degradation cost because of the division by the cost when using the true values of the parameters. The cost defined in Equation (4.12) emphasizes the stationary value of the output signal, i.e. the value which the output signal converges to eventually.

4.3 Optimal Input Signal Design: MPC of FIR system

In the examples in Section5.2an optimal input signal is designed for an FIR system controlled by an MPC controller.

Figure 4.2. MPC implemented on FIR system.

The complete system which is simulated is displayed in Figure 4.2. The notation means the same as in Figure4.1.

4.3.1 FIR system

The model structure and the transfer function from u(t) to z(t) are the same as the ones described in Section4.2.1. A disturbance is added to the system according to Figure4.2. Thus, Equation (4.5) is rewritten as

(36)

4.3.2 MPC State Space Model

The state space description of the dynamics described in Equation (4.13) is

xk+1 = Axk+ Buk,

yk = Cxk+ dk. (4.14)

The states are xk= [uk−1. . . uk−nb]

T_{. The matrices are defined as}

A=         0 0 . . . 0 0 1 0 . . . 0 0 0 1 . . . 0 0 ... ... ... ... ... 0 . . . 1 0         , B =       1 0 ... 0       , and C = [b1 b2. . . bnb],

where A is a nb-by-nb matrix and B and C are vectors of length nb. The disturbance

from time step k − 1 to time step k is obtained by

dk= yk− Cxk. (4.15)

The value of the disturbance needs to be estimated for k = 1, . . . , T , where T is the prediction horizon. The MPC algorithm assumes that all disturbances will be the same as in the initial time step, meaning that dk = d0 for k = 1, . . . , T .

4.3.3 Application Optimization Problem

The MPC controller finds the optimal sequence of input signals, {u(t)}, which minimizes the value of a particular cost function. The first input signal in the sequence is implemented on the FIR system. This procedure is then repeated for each time step k. For more details about the MPC algorithm see Section (3.2.2). The optimization problem which minimizes the cost function in MPC of FIR system can be formulated as

minimize ||Cx_MPC+ d_MPC− r_MPC||₂,

subject to xk+1 = Axk+ Buk, k= 0, . . . , T − 1, (4.16)

where || · ||2 denotes the euclidean norm. In the problem above we have xM P C =

[x0 x1. . . xT], dM P C = [d0 d1. . . dT] and rM P C = [r0 r1. . . rT]. The initial state in

each optimization problem is denoted x0. The objective of Optimization Problem

(4.16) is to minimize the effects of the disturbances and to make the measurable output, y(t), follow a reference signal, r(t). The states of the system are then updated according to

(37)

4.3. OPTIMAL INPUT SIGNAL DESIGN: MPC OF FIR SYSTEM xk+1 =         0 0 . . . 0 0 1 0 . . . 0 0 0 1 . . . 0 0 ... ... ... ... ... 0 . . . 1 0         xk+         1 0 0 ... 0         uk. (4.17)

The performance degradation cost Vapp(θ) is defined as

Vapp(θ) = 1 N N X t=1 (y(θ, t) − y(θ0, t))2. (4.18)

The different measured output signals, y(θ) and y(θ0), are obtained by simulating

the model for both θ and θ0 respectively.

4.3.5 System Identification Cost

The prediction-error-method is used to obtain the system identification matrix, see Section3.1.4, i.e. VSI(θ) = lim N →∞ 1 N N X t=1 (y(t) − ϕT_(t)θ)2_{− λ} 0, (4.19)

where N is the number of measurements used. The hessian of the system identifi-cation matrix is V_SI00(θ) = lim N →∞ 2 N N X t=1 ϕ(t)ϕT(t). (4.20)

The hessian has the simple block matrix structure

V_SI00(θ) = 2       Ru(0) Ru(1) . . . Ru(nb−1) Ru(−1) Ru(0) . . . Ru(nb−2) ... ... ... ... Ru(−(nb−1)) Ru(−(nb−2)) . . . Ru(0)       , (4.21)

where Ru(k) is the covariance function of {u(t)}, defined as

Ru(k) = lim N →∞ 1 N N X t=1 u(t)uT(t − k), (4.22) see pages 23-24 in [10].

(38)

4.3.6 System Identification Optimization Problem

It is desired to obtain the V00

SI(θ) which gives the smallest possible value on the

variance of the input signal, {u(t)}, and still ensures that the estimated system fulfills the application criterion, Vapp(θ) ≤ 1/γ. This is an optimization problem

which can be stated as

minimize Ru(0), subject to V_SI00(θ0) − γλ0χ2α(d) N V 00 app(θ0) 0, V_SI00(θ0) 0. (4.23) Note that V00

SI(θ) is calculated by using Equation (4.21). This optimization problem

is easily solved using, for instance, the mathematical programming tool MATLAB.

4.4 Optimal Input Signal Design: MPC of ARX system

In the examples described in Section5.3the optimal input signal is designed for an ARX system controlled by an MPC controller. The FIR system is a special case of the ARX formulation, see Section 3.3.2.

Figure 4.3. MPC implemented on ARX system.

The complete system which is simulated is displayed in Figure 4.3. The notation used is the same as in Figure 4.1.

4.4.1 ARX System

(39)

4.4. OPTIMAL INPUT SIGNAL DESIGN: MPC OF ARX SYSTEM

y(t) + a1y(t − 1) + · · · + anay(t − na) = b1u(t − 1) + · · ·

+bnbu(t − nb) + e(t). (4.24)

The adjustable parameters are

θ= [a1. . . ana b1. . . bnb]

T _(4.25)

and the regression vector is

ϕ(t) = [−y(t − 1) . . . − y(t − na) u(t − 1) . . . u(t − nb)]T. (4.26)

A disturbance is added to the system according to Figure 4.3. Thus, Equation (4.24) is rewritten as

y(t) + a1y(t − 1) + · · · + anay(t − na) = b1u(t − 1) + · · ·

+bnbu(t − nb) + d(t) + a1d(t − 1) + · · · + anad(t − na). (4.27)

4.4.2 MPC State Space Model

The state space description of the dynamics described in Equation (4.27) is

xk+1 = Axk+ B1uk+ B2dk,

yk = Cxk. (4.28)

The states are xk = [yk−1. . . yk−na uk−1. . . uk−nb dk−1. . . dk−na]

T_{. The matrices} are defined as A=                            

−a₁ . . . −a_n_a−1 −ana b1 . . . bnb−1 bnb 1 + a1 a2 . . . ana−1 ana

1 . . . 0 0 0 . . . 0 0 0 0 . . . 0 0 ... ... ... ... ... ... ... ... ... ... ... ... ... 0 . . . 1 0 0 . . . 0 0 0 0 . . . 0 0 0 . . . 0 0 0 . . . 0 0 0 0 . . . 0 0 0 . . . 0 0 1 . . . 0 0 0 0 . . . 0 0 ... ... ... ... ... ... ... ... ... ... ... ... ... 0 . . . 0 0 0 . . . 1 0 0 0 . . . 0 0 0 . . . 0 0 0 . . . 0 0 0 0 . . . 0 0 0 . . . 0 0 0 . . . 0 0 1 0 . . . 0 0 0 . . . 0 0 0 . . . 0 0 0 1 . . . 0 0 ... ... ... ... ... ... ... ... ... ... ... ... ... 0 . . . 0 0 0 . . . 0 0 0 0 . . . 1 0                             ,

(40)

B1 = [0 0 . . . 0 1 0 . . . 0 0 0 0 . . . 0]T,

B2 = [0 0 . . . 0 0 0 . . . 0 1 0 0 . . . 0]T

and

C = [−a1 − a2. . . − ana b1 b2. . . bnb 1 + a1 a2. . . ana] (4.29)

where A is a n-by-n matrix and B1, B2 and C are vectors of length n with n =

2na+nb. The nonzero elements of B1and B2occur at row na+1 and row na+nb+1.

The disturbance from time step k − 1 to time step k is obtained by

dk= yk− Cxk+ dk−1. (4.30)

The value of the disturbance one time step ahead needs to be estimated. The MPC algorithm assumes that the next disturbance will be the same as the previous one, meaning that dk= dk−1. This explains the appearance of the element (1+a1) which

may look out of place in the matrix A and the vector C.

4.4.3 Application Optimization Problem

The optimization problem which minimizes the cost function in the MPC of ARX system can be formulated as

minimize ||CxMPC− rMPC||2,

subject to xk+1 = Axk+ B1uk+ B2dk, k= 0, . . . , T − 1, (4.31)

where the same notation as in Optimization Problem (4.16) is used. The objective of Optimization Problem (4.31) is to minimize the effects of the disturbances and to make the measurable output, y(t), follow a reference signal, r(t). The states of the system are then updated according to

xk+1 =                           0 . . . 0 0 1 . . . 0 0 ... ... ... ... 0 . . . 1 0 0 . . . 0 0 1 . . . 0 0 ... ... ... ... 0 . . . 1 0 0 . . . 0 0 1 . . . 0 0 ... ... ... ... 0 . . . 1 0                           xk+                           1 0 ... 0 0 0 ... 0 0 0 ... 0                           yk+                           0 0 ... 0 1 0 ... 0 0 0 ... 0                           uk+                           0 0 ... 0 0 0 ... 0 1 0 ... 0                           dk. (4.32)

(41)

4.4. OPTIMAL INPUT SIGNAL DESIGN: MPC OF ARX SYSTEM

The performance degradation cost Vapp(θ) is defined in the same way as for the FIR

system, see Equation (4.18).

4.4.5 System Identification Cost

The prediction-error-method is used as in Section4.3.5, i.e. the hessian of the system identification cost is V_SI00(θ) = lim N →∞ 2 N N X t=0 ϕ(t)ϕT(t). (4.33)

The hessian has the simple block matrix structure

V_SI00(θ) = 2 Ry −Ryu −R_uy R_u ! , (4.34) with Ry =       Ry(0) Ry(1) . . . Ry(nb−1) Ry(−1) Ry(0) . . . Ry(nb−2) ... ... ... ... Ry(−(nb−1)) Ry(−(nb−2)) . . . Ry(0)       , (4.35) R_u =       Ru(0) Ru(1) . . . Ru(nb−1) Ru(−1) Ru(0) . . . Ru(nb−2) ... ... ... ... Ru(−(nb−1)) Ru(−(nb−2)) . . . Ru(0)       , (4.36) and R_uy=       Ruy(0) Ruy(1) . . . Ruy(nb−1) Ruy(−1) Ruy(0) . . . Ruy(nb−2) ... ... ... ... Ruy(−(nb−1)) Ruy(−(nb−2)) . . . Ruy(0)       , (4.37)

where Ry(k) is the covariance function of {y(t)}, Ru(k) is the covariance function of

{u(t)}, Ruy(k) and Ryu(k) are the covariance functions of {y(t)} and {u(t)}. The

(42)

Ruy(k) = lim N →∞ 1 N N X t=1 u(t)yT(t − k), (4.38)

The matrix structure of V00

SI(θ) for the FIR system only consists of the lower right

corner, Ru, see Section 4.3.5.

4.4.6 System Identification Optimization Problem

As for the FIR system, it is desired to obtain the V00

SI(θ) which gives the smallest

possible value on the variance of the input signal, {u(t)}, and still ensures that the estimated system fulfills the application criterion, Vapp(θ) ≤ 1/γ. This is an

optimization problem which can be stated as minimize Ru(0), subject to V_SI00(θ0) − γλ0χ2α(d) N V 00 app(θ0) 0, V_SI00(θ0) 0, (4.39) Note that V00

SI(θ0) is calculated by using Equation (4.34). This optimization problem

is easily solved if na= 0, i.e. if the system has a FIR structure. If not, the problem

has to be translated to an equivalent one in the frequency domain.

4.4.7 Frequency Domain

The spectral density is denoted by Φ(ω). It is defined as Φ(ω) =

∞

X

k=−∞

R(k)e−ikω, (4.40)

where R(k) is the covariance function, which is defined as

R(k) = 1

2π Z π

−πΦ(ω)e

ikω_dω. _(4.41)

Under certain conditions, see pages 40-41 in [10], one can show that

Φy(ω) = |G0(eiω)|2Φu(ω) + λ0 (4.42)

and

Optimal Input Signal Design and MPC of Nonlinear Dynamical Systems