Mikles J., Fikar M. Process Modelling, Identification, and Control

Подождите немного. Документ загружается.

422 9 Predictive Control

y(k + N

+ i)=w(k + N

),i=1,...,n (9.80)

and n is the dimension of the system state vector.

Although the output constraints are added to the original formulation, the

solution in the unconstrained case can still be found analytically.

Theorem 9.6 (SIORHC). Let the system polynomials ΔA and B be co-

prime. If λ>0 and deg(B) ≤ deg(A)+1 then provided that

≥ n = max(deg(A)+1, deg(B))

then

• the SIORHC control law is unique;

• SIORHC stabilises the plant, and, irrespective of deg(A), deg(B),λ,for

= n

yields a state dead-beat closed-loop system;

• whenever stabilising, SIORHC yields asymptotic rejection of constant dis-

turbances and oﬀset free closed-loop system.

Derivation

As usual in GPC, consider the predicted output to be of the form

ˆy = G˜u + y

(9.81)

where all vectors are stacked from k +1 to k + N

. After this time, up to

k + N

+ n, output predictions are constrained to be equal to the setpoint and

future control increments to zero

¯w =

G˜u + ¯y

(9.82)

where the vector ¯y

denotes the free response between k+N

+1 and k+N

and ¯w is vector of w(k + N

) of the corresponding dimension.

The cost function is as usual of the form

I =(ˆy − w)

(ˆy − w)+λ˜u

˜u (9.83)

subject to the constraint (9.82). This can be solved analytically by the La-

grange multipliers.

Let us denote

h = ¯w − ¯y

and h = w − y

. Hence

I =(G˜u − h)

(G˜u −h)+λ˜u

˜u + x

(

G˜u −

h) (9.84)

Partial derivatives of I with respect to ˜u, x (x are the Lagrange multipliers)

are zero. We obtain a system linear equations as follows

9.5 Stability Results 423



2(G

G + λI)

G 0



˜u







(9.85)

The block matrix inversion formula states



−1



−1



A + ADΔCA −ADΔ

−ΔCA −Δ



,Δ

−1

= B − CAD (9.86)

Therefore, the future control increment vector is given as

˜u =

G(I −

G)G

h +

h (9.87)

where

G =



G + λI



−1

(9.88)

Q =





−1

(9.89)

SGPC

Another method that can be shown to be equivalent to the preceding methods

is SGPC. Its advantages include more eﬃcient computational implementation

and better numerical robustness.

In this approach GPC is invoked after application of the stabilising feed-

back control law

Y (q

−1

)Δu(t)=c(t) −X(q

−1

)y(t) (9.90)

where the polynomials X, Y are calculated as the dead-beat controller from

the Diophantine equation

ΔAY + BX = 1 (9.91)

and where c(t) denoted the reference signal for the closed-loop system that

forms the vector of optimised variables. The dead-beat controller results in

the control and output predictions of the form

y(t)=B(q

−1

)c(t) (9.92)

u(t)=A(q

−1

)c(t) (9.93)

These can be simulated forward in time to give the vectors of future output

and control predictions and treated in the same way as in GPC.

Hence, this methods optimises future reference trajectory subject to (ter-

minal) constraint that this trajectory should be equal to the desired setpoint

after some horizon.

It can be shown that SGPC is equivalent to the pole-placement method

with the controller polynomials Y

given from the Diophantine equation

ΔAY

+ BX

= P

(9.94)

and its stability depends on the roots of the P

polynomial.

Theorem 9.7 (SGPC). For N

≥ deg(A)+1+N

,whereN

is the number

of reference points optimised, is SGPC stable for any N

424 9 Predictive Control

YKPC

The last approach within ﬁnite horizon formulation is the predictive control

algorithm based on the Youla Kuˇcera parametrisation of all stabilising con-

trollers YKPC.

As in the previous approach, in the ﬁrst step a nominal controller with

pole-placement technique is calculated and a nominal controller is given as a

solution of 2 Diophantine equations

ΔA + Q

B = M (9.95)

ΔS

+ BR

= M (9.96)

where M is the desired closed-loop polynomial and the two degree of free-

dom controller with integral action is deﬁned as Q/ΔP (feedback) and R/ΔP

(feedforward).

The minimum degree controller P

, Q

, R

only serves as a basis to ﬁnd

an expression for the set of all stabilizing controllers. Among these controllers

the one is chosen that minimises the GPC cost function.

The expression of such controllers (Youla-Kuˇcera parametrisation) is as

follows:

Theorem 9.8 (YK Controllers). A controller (P (z),Q(z),R(z)) gives rise

to the closed-loop denominator matrix M (z) if and only if it can be

expressed as

P = P

+ ZB (9.97)

Q = Q

− ZAΔ (9.98)

R = R

+ ΔX (9.99)

X, Z are assumed to be polynomials for simplicity. Their coeﬃcients form

the vector of the optimised parameters.

Stability is proved as in the previous approaches via terminal constraint.

It is interesting to note, that in this approach the state terminal constraint

does not have to be speciﬁed and is implicitly assured.

Theorem 9.9 (Choice of horizons). Let n = max(deg(X), deg(Z)), N

and let the horizons be equal or greater than

= deg(B)+n (9.100)

= deg(AΔ)+n (9.101)

Further assume that the sequences w, d (reference, disturbance) are bounded.

Then the unconstrained YKPC is uniformly asymptotically stable.

This controller is time-varying in spite of the fact that the system is as-

sumed to be time-invariant.

9.5 Stability Results 425

9.5.3 Inﬁnite Horizons

Another line of research has been focused on reformulation of the basic GPC

method when N

are inﬁnity. Of course, if such a method can be imple-

mented, stability problems disappear. However, a number of the optimised

parameters (future control moves) is also inﬁnity and the original problem

is untractable. Therefore, several suboptimal algorithms have emerged. The

basic principle of all of them is to leave N

= ∞ but to play with N

or with

its equivalents.

Rawlings a Muske developed a method in state-space formulation where the

number of control moves N

is ﬁnite. The feedback gain is calculated via the

recursive ARE.

Theorem 9.10 (Stable plants). For stable system matrix A and N

≥ 1 is

the receding horizon controller stabilising.

Theorem 9.11 (Unstable plants). For stabilisable plant (A, B) with r un-

stable modes and N

≥ r is the receding horizon controller stabilising.

Constrained control can also be dealt with in this approach. The require-

ment added to the previous theorems is that the initial state at time k is

feasible (within constraints).

The SGPC and YKPC methods can be modiﬁed to use both input and output

horizons inﬁnite. The SGPC approach utilises the ﬁnite reference sequence as

the vector of optimised variables. The solution is found via Lyapunov equation.

The YKPC method utilises coeﬃcients of the Youla-Kuˇcera polynomials

as the optimised variables. It is shown that in the unconstrained case the

optimal predictive controller coincides with the nominal pole-placement con-

troller whose poles are calculated via spectral factorisation equation – hence

it is the standard LQ controller. If the constraints are active, piece-wise linear

controller results.

9.5.4 Finite Terminal Penalty

The third approach to MBPC stability is to adopt a ﬁnite input and state

horizon with a ﬁnite terminal weighting matrix. This is equivalent to the

condition that the terminal state has to be within some neighbourhood of

the origin. Compared to the previous approaches when it had to be exactly

zero, here the state can be in such a neighbourhood of the origin that is

asymptotically stable.

With a state-space formulation

x(k +1)=A(k)x(k)+B(k)u(k) (9.102)

is the cost given as

426 9 Predictive Control

I = x

(k + N)P (k + N)x(k + N)+



i=0



(k + i)Q(k + i)x(k + i)

+ u

(k + i)R(k + i)u(k + i)



(9.103)

If P (k + N) →∞then the terminal state constraint approach results. How-

ever, also some “smaller” terminal penalty matrix P can give a stabilising

receding horizon control. Equivalently, in terms of state constraints, we will

speak about stabilising state terminal sets.

Theorem 9.12. Assume that the terminal weighting matrix P (k+1) satisﬁes

the following matrix diﬀerence inequality for some matrix H(k)

P (k) ≥ F

(k)P (k +1)F (k)+Q(k)+H

(k)R(k +1)H(k), ∀k ∈ [N,∞)

(9.104)

where

F (k)=A(k)+B(k)H(k) (9.105)

Further suppose that Q(k) is positive deﬁnite. Then the receding horizon con-

trol which stems from the optimisation problem minimising the performance

index I asymptotically stabilises the system. In addition, if A(k), Q(k),and

P (k) are bounded above ∀k ≥ 0, then the receding horizon control exponen-

tially stabilises the system.

9.6 Explicit Predictive Control

In receding horizon control, a ﬁnite-time optimal control problem is solved at

each time step to obtain the optimal input sequence. Subsequently, only the

ﬁrst element of that sequence is applied to the system. At the next time step,

the state is measured and the procedure is repeated from the beginning. The

input sequence can be computed by solving an optimisation problem on-line at

each time step. Hence the complexity of the underlying optimisation problem

limits the minimal admissible sampling time of the plant.

Alternatively, it is possible to solve the optimisation problem oﬀ-line as

a multi-parametric program. Then, the on-line eﬀort reduces to ﬁnding the

correct feedback law entry in a lookup table, allowing MPC to be applied to

systems with very low sampling times.

9.6.1 Quadratic Programming Deﬁnition

Consider optimal control problems for a discrete-time linear, time-invariant

(LTI) system

9.6 Explicit Predictive Control 427

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

with A[n ×n]andB[m ×n].

Now consider the constrained ﬁnite-time optimal control problem

min I(x(k)) = x

(k + N)P (k + N)x(k + N)



i=0



(k + i)Q(k + i)x(k + i)

+ u

(k + i)R(k + i)u(k + i)



(9.107a)

subj. to x(k + i) ∈ X, u(k + i) ∈ U,i∈{0,...,N − 1}, (9.107b)

x(k + N) ∈T

set

, (9.107c)

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

(9.107d)

where we assume that states and control are subject to constraints (9.107b).

The terminal set constraint (9.107c) is an additional constraint which is often

added to obtain certain properties (i.e. stability and constraint satisfaction).

As future state predictions are constrained by (9.106), we can recursively

substitute for them yielding

x(k + i)=A

x(k)+

i−1



j=0

Bu(k + i − 1 − j) (9.108)

Thus, optimal solution of problem (9.107) can be reformulated as

∗

)=x

+ min



+ x



(9.109a)

subj. to GU

≤ W + Ex

(9.109b)

where the column vector U

=[u

,...,u

N−1

]

is the optimisation vector

and H, F , Y , G, W , E can easily be obtained from the original formulation.

The reformulated problem (9.109) is a standard quadratic programming

formulation that was also obtained in previous sections. Taking any initial

value x

, optimal future control trajectory U

(x) can be found from which

only the ﬁrst element is used to close the loop.

9.6.2 Explicit Solution

In explicit solution of (9.109) we use the so-called multi-parametric program-

ming approach to optimisation. In multi-parametric programming, the objec-

tive is to obtain the optimiser U

for a whole range of parameters x

, i.e.

to obtain U

(x) as an explicit function of the parameter x.Thetermmulti

is used to emphasise that the parameter x (in our case the actual state vec-

tor x

) is a vector and not a scalar. If the objective function is quadratic

428 9 Predictive Control

in the optimisation variable U

, the terminology multi-parametric Quadratic

Program (mp-QP) is used.

In this formulation, it is useful to deﬁne

z = U

+ H

−1

(9.110)

and to transform the formulation to problem, where the state vector x

apears

only in constraints

∗

) = min





(9.111a)

subj. to Gz ≤ W + Sx

(9.111b)

where S = E + GH

−1

An mp-QP computation scheme consists of the following three steps:

1. Active Constraint Identiﬁcation: A feasible parameter

x is determined and

the associated QP (9.111) is solved. This will yield the optimiser z

and active constraints A(

x) deﬁned as inequalities that are active at

solution, i.e.

x)={i ∈J |G

(i)

z = W

(i)

+ S

(i)

x}, J = {1, 2,...,q}, (9.112)

where G

(i)

, W

(i)

,andS

(i)

denote the i-th row of the matrices G, W ,

and S, respectively, and q denotes the number of constraints. The rows

indexed by the active constraints A(

x) are extracted from the constraint

matrices G, W and S in (9.111) to form the matrices G

, W

and S

2. Region Computation: Next, it is possible to use the Karush-Kuhn-Tucker

(KKT) conditions to obtain an explicit representation of the optimiser

(x) which is valid in some neighborhood of

x. These are for our prob-

lem deﬁned as

Hz + G

λ = 0 (9.113a)

(Gz − W − S

x)=0 (9.113b)

λ ≥ 0 (9.113c)

Gz ≤ W + S

x (9.113d)

Optimised variable z can be solved from (9.113a)

z = −H

−1

λ (9.114)

Condition (9.113b) can be separated into active and inactive constraints.

For inactive constraints holds λ

= 0. For active constraints are the cor-

responding Lagrange multipliers λ

positive and inequality constraints

are changed to equalities. Substituting for z from (9.114) into equality

constraints gives

9.6 Explicit Predictive Control 429

−G

−1

+ W

+ S

x = 0 (9.115)

and yields expressions for active Lagrange multipliers

= −(G

−1

)

−1

+ S

x) (9.116)

The optimal value of optimiser z and optimal control trajectory U

are

thus given as aﬃne functions of

z = −H

−1

)

−1

+ S

x) (9.117)

= z − H

−1

= −H

−1

)

−1

+ S

x) −H

−1

= F

x + G

(9.118)

where

= H

−1

)

−1

− H

−1

(9.119)

= H

−1

)

−1

(9.120)

In a next step, the set of states is determined where the optimiser U

(x)

satisﬁes the the same active constraints and is optimal. Such a region is

characterised by two inequalities (9.113c), (9.113d) and is written com-

pactly as H

x ≤ K

where



G(F

+ H

−1

) −S

−1

)

−1



(9.121)



W − GG

−(G

−1

)

−1



(9.122)

3. State Space Exploration: Once the controller region is computed, the algo-

rithm proceeds iteratively in neighbouring regions until the entire feasible

state space is covered with controller regions.

After completing the algorithm, the explicit model predictive controller

consists of deﬁnitions of multiple state regions with diﬀerent aﬃne control

laws. Its actual implementation reduces to search for an active region of states

and calculation of the corresponding control.

9.6.3 Multi-Parametric Toolbox

The Multi-Parametric Toolbox (MPT) is a free, open-source, Matlab-based

toolbox for design, analysis and deployment of optimal controllers for con-

strained linear, nonlinear and hybrid systems. The toolbox oﬀers a broad spec-

trum of algorithms developed by the international community and compiled

in a user friendly and accessible format: starting from diﬀerent performance

objectives (linear, quadratic, minimum time) to the handling of systems with

430 9 Predictive Control

persistent additive and polytopic uncertainties. Users can add custom con-

straints, such as polytopic, contraction, or collision avoidance constraints, or

create custom objective functions. Resulting optimal control laws can either

be embedded into applications in the form of C code, or deployed to target

platforms using the Real Time Workshop.

Although the primal focus of MPT is on design of optimal controllers based

on multi-parametric programming, it can be used to formulate and solve on-

line MPC problems as well. The toolbox contains eﬃcient algorithms to solve

linear and quadratic problems parametrically and features a rich library for

performing operations on convex polytopes and non-convex unions thereof.

Example 9.2:

www

To illustrate the capabilities of the Multi-Parametric Toolbox, we consider

a double integrator whose transfer function representation is given by

P (s)=

With a sampling time T

= 1 second, the corresponding state-space form

can be written as

x(k +1)=





x(k)+



0.5



u(k)

y(k)=x

(k)

We want to design an explicit optimal state-feedback law which minimizes

the performance index (9.107a) with N =5,P = 0, Q = I,andR =1.

System states and the control signals are subject to constraints x(k) ∈

[−1, 1] × [−1, 1] and u(k) ∈ [−1, 1], respectively. In order to solve the

problem (9.107) explicitly, one has ﬁrst to describe the dynamical model

of the plant:

model.A = [1 1; 0 1];

model.B = [1; 0.5];

model.C = [0 1];

model.D = 0;

along with the system constraints:

model.umin = -1;

model.umax = 1;

model.xmin = [-1; -1];

model.xmax = [1; 1];

Next, parameters of the performance index have to be speciﬁed:

cost.N = 5;

cost.Q = [1 0; 0 1];

cost.R = 1;

cost.P_N = 0;

cost.norm = 2;

9.7 Tuning 431

Finally, the explicit optimal state-feedback control law can be calculated

by executing

controller = mpt_control(model, cost)

The obtained explicit control law is deﬁned by

u =

⎧

⎪

⎨

⎪

⎩



−0.52 −0.94



x, if

⎛

⎜

⎝

0.48 0.88

−0.48 −0.88

−10

0 −1

⎞

⎟

⎠

x ≤

⎛

⎜

⎝

0.93

⎞

⎟

⎠

(Region #1)

−1, if

⎛

⎝

−0.48 −

0.88

⎞

⎠

x ≤

⎛

⎝

−0.93

⎞

⎠

(Region #2)

1, if

⎛

⎝

−10

0 −1

0.48 0.88

⎞

⎠

x ≤

⎛

⎝

−0.93

⎞

⎠

(Region #3)

and it can be plotted using the command

plot(controller)

which will show the regions of the state-space over which the optimal

control law is deﬁned, as illustrated in Figure 9.3.

9.7 Tuning

Let us recall the GPC parameters: horizons N

, control weighting λ,

and output weighting polynomial P . The eﬀects of a change of the parame-

ters are strongly coupled and the strategies dealing with adjustment of GPC

parameters usually adjust only one parameter while all others are at some

default values.

9.7.1 Tuning based on the First Order Model

The tuning strategy based on the analysis of the ﬁrst order system with time

delay is given as follows: The strategy is given as follows:

1. Approximate the process dynamics with the ﬁrst order plus dead time

model:

G(s)=

Ts+1

−T

(9.123)