Stable multivariable model based predictive control

Sharon A. Heise^{∗} and Jan M. Maciejowski Department of Engineering Cambridge University Cambridge CB2 1PZ England

Fax: +44 223 332662

                                 email: sah@eng.cam.ac.uk
                                    jmm@eng.cam.ac.uk

33rd CDC

Abstract

The stability of model based predictive control (MBPC) is studied in a multiple-input, multiple-output framework in the presence of active state and input constraints, and implemented in both an input-output and a state space format. The constrained generalized predictive control method of Rossiter and Kouvaritakis is extended to the multivariable case. In this approach, stability is established by demonstrating the existence of a monotonically decreasing cost function.

Keywords: Stability; Constraints; Predictive Control.

1 Introduction

Stability in the presence of constraints has long been an issue of concern in model based predictive control (MBPC), largely due to the nonlinear nature of the overall control scheme. Analysis of system stability is straightforward and may be considered in terms of a contraction mapping condition [16] or, given the piecewise linear time invariant characteristic of MBPC in which an LTI controller may be obtained when a particular set of constraints is active, a small gain related condition [6]. These analysis procedures are inherently conservative, however, so it is desirable to obtain a design procedure which guarantees stability.

Within the last few years, it has been found that treating the MBPC problem on an infinite horizon results in a closed loop stable system in the presence of constraints as long as a feasible solution to the optimization problem exists, at least for state feedback [10], [11]. At the same time, it has often been assumed that there are only a few highly restrictive stability results to be obtained from a finite horizon formulation. The method of stable generalized predictive control (SGPC) [7], however, provides a novel means of developing control strategies with guaranteed stability for singleinput, single-output systems on a finite horizon. SGPC involves the incorporation of an internally stabilizing feedback loop within an optimization and prediction scheme, resulting in simple closed loop relations which yield finite closed loop impulse responses, thus ensuring the resulting optimized cost is monotonically decreasing. This monotonically decreasing cost property is then used to establish closed loop stability. Furthermore, the method of constrained stable generalized predictive control (CSGPC) [13] extends SGPC to address constraints by augmenting the SGPC cost and applying a mixed weights least squares optimization. The resulting design algorithm achieves optimality, subject to constraint satisfaction, while still retaining a monotonically decreasing cost property for single-input, single-output systems in input-output form.

∗This author for correspondence.

The purpose of this paper is to extend the CSGPC method to address multiple-input, multiple output systems, first in input-output and then in state space form. In the next section we will establish the multivariable equivalent of the stabilizing inner feedback loop, show that the monotonically decreasing cost property holds in the multivariable case for specified prediction, control and reference horizons, and illustrate the implementation of the optimization strategy. The state space implementation of such a method will then be shown. Feasibility and convergence issues will not be addressed, and we will concentrate on stability.

2 MIMO stabilizing inner feedback loop

In this section we follow the approach of [13], generalizing as necessary to the multivariable case and removing the assumption that the control horizon is the same as the prediction horizon. A plant G(z) ∈ RL_{∞} is assumed to have a double coprime factorization:

G(z) = N_{G}(z)M_{G}(z)^{−1} = M̃_{G}(z)^{−1}Ñ_{G}(z)  (1)

where M_{G}(z), N_{G}(z) ∈ RH_{∞}. The assumption of double coprimeness implies that there exist matrices transfer function X, Y, X̃ and Ỹ ∈ RH_{∞} such that

[ \begin{matrix}Ỹ & X̃ \\ −Ñ_{G} & M̃_{G} \end{matrix} ] [ \begin{matrix}M_{G} & −X \\ N_{G} & Y \end{matrix} ] = I  (2)

Then the solution of the Bezout identity

X̃N_{G} + Ỹ M_{G} = I  (3)

yields the set of all proper controllers achieving internal stability for G [5], which may be parameterized by

K = −(Ỹ − QÑ_{G})^{−1}(X̃ + QM̃_{G})  (4)

for Q ∈ RH_{∞}. It will be convenient to let M_{G} and N_{G} take the form:

M_{G}(z) = [ \begin{matrix}m_{11}(z)Δ & … & m_{1ξ}(z)Δ \\ ⋮ & ... & ⋮ \\ m_{ξ1}(z)Δ & … & m_{ξξ}(z)Δ \end{matrix} ] = M(z)  (5)

N_{G}(z) = [ \begin{matrix}z^{−1}n_{11}(z) & … & z^{−1}n_{1ξ}(z) \\ ⋮ & ... & ⋮ \\ z^{−1}n_{η1}(z) & … & z^{−1}n_{ηξ}(z) \end{matrix} ] = N(z)  (6)

where Δ is the differencing operator

Δ = 1 − z^{−1}

and each element of N(z) and M(z) is finite impulse response (FIR) filter of the form

n_{ij}(z) = n_{ij}^{0} + n_{ij}^{1}z^{−1} + … + n_{ij}^{n}z^{−n}  (7)
m_{ij}(z) = m_{ij}^{0} + m_{ij}^{1}z^{−1} + … + m_{ij}^{n}z^{−n}

It will be helpful in the sequel to note that M(z) and N(z) also take the form

M(z) = T_{u}^{0} + T_{u}^{1}z^{−1} + … + T_{u}^{n+1}z^{−(n+1)}
N(z) = T_{y}^{0}z^{−1} + T_{y}^{1}z^{−2} + … + T_{y}^{n}z^{−(n+1)}  (8)

where T_{y}^{τ} ∈ ^{η×ξ} is a matrix with ijth element n_{ij}^{τ} and T_{u}^{τ} ∈ ^{ξ×ξ} consists of ijth elements m_{ij}^{τ} − m_{ij}^{(τ−1)}, with the special cases of T_{u,ij}^{0} = m_{ij}^{0} and T_{u,ij}^{n} = −m_{ij}^{n}.

No special assumptions are made about Ñ_{G} or M̃_{G} other than those already stated, so long as (2) is satisfied. The Bezout identity (3) to be solved then becomes:

X̃N + Ỹ M = I  (9)

The stabilizing feedback loop is then realized as in figure 1, where v(k) ∈ ^{ξ} is an artificial reference signal, y(k) ∈ ^{η} is the system output, and Δu(k) ∈ ^{ξ} is the control input.

[figure: y; Δ; u; K; G; v; M^{−1}; N; (Ỹ − QÑ_{G})^{−1}; X̃ + QM̃_{G}; −]

Figure 1: Stabilizing inner feedback loop

From the above figure and equations (4) and (9), we find that the following closed loop relations hold:

y(z) = N(z)v(z)  (10)
Δu(z) = M(z)v(z)  (11)

This formulation implies that the plant operates on the control increments Δu, rather than the actuator inputs u, which introduces integral action into the feedback loop. Note that N(z) and M(z) are FIR operators.

2.1 Prediction

The future system behavior may be simulated forward in time across a prediction horizon p and a control horizon q by first expanding (10) and (11) in terms of the difference equations implied by (5) and (6):

y(k + l) = ∑_{τ=0}^{n} T_{y}^{τ} v(k + l − τ − 1)  l = 1 … p  (12)

Δu(k + l) = ∑_{τ=0}^{n+1} T_{u}^{τ} v(k + l − τ)  l = 0 … q − 1  (13)

Define the vector of future system outputs y and future input increments Δu as

ỹ(k) = [ y(k + 1)^{T} … y(k + p) ]^{T}  (14)
Δũ(k) = [ Δu(k)^{T} … Δ(k + q − 1) ]^{T}  (15)

where Δu(k + l) = 0 for l ≥ q. Also define r as an optimization horizon, which is not necessarily the same as the control horizon as assumed in [13], and denote the past and future values of the reference signal v(k) by

ṽ(k) = [ v(k)^{T} … v(k + r − 1)^{T} ]^{T}  (16)
ṽ^{b}(k) = [ v(k − 1)^{T} … v(k − n)^{T} ]^{T}  (17)
ṽ^{∞}(k) = [ v^{∞}(k + r)^{T} … v^{∞}(k + p − 1) ]^{T}  (18)

In general we assume that p ≥ q ≥ r. As will be seen in the next section, a convenient choice for the far future reference signal v^{∞}(k) is v^{∞}(k) = v^{∞}, where v^{∞} is the solution of the system of equations

[ \begin{matrix}n_{11}(1) & … & n_{1ξ}(1) \\ ⋮ & ... & ⋮ \\ n_{η1}(k) & … & n_{ηξ}(1) \end{matrix} ] [ v_{1}^{∞} ; ⋮ ; v_{ξ}^{∞} ] = [ s_{1}^{∞} ; ⋮ ; s_{η}^{∞} ]  (19)

where s^{∞} is the steady state value of a desired output setpoint trajectory. Then the future system behavior may be described in terms of prediction equations:

ỹ(k) = Γ_{y}ṽ(k) + Ξ_{y}ṽ^{b}(k) + X_{y}ṽ^{∞}(k)  (20)
Δũ(k) = Γ_{u}ṽ(k) + Ξ_{u}ṽ^{b}(k) + X_{u}ṽ^{∞}(k)  (21)

where Γ_{y} ∈ ^{pη×rξ} and X_{y} ∈ ^{pη×(p−r)ξ} are Toeplitz matrices in T_{y}^{τ}, Γ_{u} ∈ ^{qξ×rξ} and X_{u} ∈ qξ×(p−r)ξ are Toeplitz matrices in T_{u}^{τ}, Ξ_{y} ∈ ^{pη×nξ} is a Hankel matrix in T_{y}^{τ} and Ξ_{u} ∈ ^{qξ×nξ} is a Hankel matrix in T_{u}^{τ}. Specifically,

Γ_{y} = [ \begin{matrix}T_{y}^{0} & 0 & … & … & … & … & 0 \\ T_{y}^{1} & T_{y}^{0} & 0 & … & … & … & 0 \\ ⋮ & ⋮ & ... & ⋮ & ⋮ & ⋮ & ⋮ \\ T_{y}^{n} & T_{y}^{n−1} & … & T_{y}^{0} & … & … & 0 \\ 0 & T_{y}^{n} & T_{y}^{n−1} & … & … & … & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & … & … & … & T_{y}^{1} & T_{y}^{0} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & … & … & … & … & … & 0 \end{matrix} ]  (22)

Ξ_{y} = [ \begin{matrix}T_{y}^{1} & T_{y}^{2} & … & T_{y}^{n} \\ T_{y}^{2} & … & T_{y}^{n} & 0 \\ ⋮ & ... & ⋮ & ⋮ \\ T_{y}^{n} & 0 & … & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & … & … & 0 \end{matrix} ]  (23)

X_{y} = [ \begin{matrix}0 & … & … & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ T_{y}^{0} & 0 & … & 0 \\ T_{y}^{1} & T_{y}^{0} & … & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & … & T_{y}^{1} & T_{y}^{0} \end{matrix} ]  (24)

The terms corresponding to the Δũ equation are similar, with the addition of a T^{n+1} term in the matrices, and so they will not be included here.

2.2 Stability without constraints

The goal of the model based predictive control problem is to minimize the deviation of the future system outputs from a specified setpoint trajectory, subject to both a penalty on control power and constraints. Control system performance in meeting this goal is qualified by the following performance index, which is to be minimized over the reference signal ṽ(k):

J(k) = ∑_{l=1}^{p} ∥y(k + l) − s(k + l)∥_{Q(l)}^{2} + ∑_{l=0}^{q} ∥Δu(k + l)∥_{R(l)}^{2}  (25)

where s(k) ∈ ^{η} is a vector of system setpoints. Q(l) and R(l) are positive definite and positive semi-definite weighting matrices, respectively, which may vary across the control horizon q and the prediction horizon p, but are assumed to be independent of k. Q(l) provides an extra degree oftuning freedom beyond that normally assumed.

We will first prove the stability of the inner feedback loop without consideration of constraints by generalizing two theorems of [7] using [17]; the first addresses the setpoint tracking error only, while the second adds a penalty on control power.

Theorem 1 Let R(l) = 0 and assume p > r+α, where α is the highest degree of any n_{ij}(z). Then, the system of figure (1) is stable for all r.

Proof Without affecting the proof, we can assume the setpoints s̃(k) are constant (step commands), and let ṽ^{∗}(k) be the optimal reference signal that minimizes (25) at time k. Given a ṽ^{∞}(k) which solves the system of equations (19), then y(k) = N(z)v(k) will reach its setpoint value within p time steps or less, as long as p > r + α. Now ṽ^{∗}(k) is a feasible solution to the optimization at time k + 1, but not necessarily optimal. If ṽ^{∗}(k) where implemented at k + 1, the cost would be

J(k + 1) = J^{∗}(k) − [y(k + 1) − s(k + 1)]^{T} Q(1) [y(k + 1) − s(k + 1)] < J^{∗}(k)  (26)

However, at time k + 1 we will reoptimize, and thus

J^{∗}(k + 1) ≤ J(k + 1) < J^{∗}(k)  (27)

In addition, feasibility of the optimization implies that J_{1} is finite. Now we have

J(k + 1) + ∑_{l=1}^{k} [y(l + 1) − s(l + 1)]^{T} Q(l) [y(l + 1) − s(l + 1)] + ≤ J_{1} < ∞  (28)

Thus the cost is finite and monotonically decreasing. In the limit, then, y will reach its target value and stability will be achieved.

□

Theorem 2 Assume R(l) ≥ 0 and p > r + α where α is the highest degree of any n_{ij}(z). Then, the system of figure (1) is stable for all r.

Proof The same proof as in theorem 1 applies, with slight modification. Now when ṽ^{∗}(k) is implemented at k + 1, the cost is

J(k + 1) = J^{∗}(k) − [y(k + 1) − s(k + 1)]^{T} Q(1) [y(k + 1) − s(k + 1)] −
Δu(k)^{T} R(0)Δu(k) < J^{∗}(k)  (29)

And the monotonically decreasing cost becomes

J(k + 1) + ∑_{l=1}^{k} {[y(l + 1) − s(l + 1)]^{T} Q(l) [y(l + 1) − s(l + 1)] +
+ Δu(l)^{T} R(l)Δu(l)} ≤ J_{1} < ∞  (30)

Thus y will reach its setpoint value in the limit and stability is established.

□

2.3 Least squares cost formulation

The performance index of equation (25) may be expanded in terms of l by using equations (20) and (21). Define

Q = diag [ Q(1) … Q(p) ]
R = diag [ R(0) … R(q) ]  (31)

And write the cost function in terms of ṽ as:

J(k) = ṽ(k)^{T} Sṽ(k) + 2 {[ỹ_{c}(k) − s̃(k)]^{T} QΓ_{y} + Δũ_{c}(k)^{T} RΓ_{u}} ṽ(k) +
[ỹ_{c}(k) − s̃(k)]^{T} Q [ỹ_{c}(k) − s̃(k)] + Δũ_{c}(k)^{T} RΔũ_{c}(k)  (32)

where

S = Γ_{y}^{T} QΓ_{y} + Γ_{u}^{T} RΓ_{u}

and ỹ_{c}(k) and Δũ_{c}(k) are the contributions to ỹ(k) and Δũ(k) that are independent of the future reference signals ṽ(k). That is, from (20) and (21),

ỹ_{c}(k) = Ξ_{y}ṽ^{b}(k) + X_{y}ṽ^{∞}(k)
Δũ_{c}(k) = Ξ_{u}ṽ^{b}(k) + X_{u}ṽ^{∞}(k)  (33)

As in [13], it is convenient to complete the squares to form a Euclidean norm on J, denoted by ∥·∥, since the addition of constants to the cost will have no effect on its minimization. Let

ṽ^{0}(k) = S^{−1} {(QΓ_{y})^{T} [s̃(k) − ỹ_{c}(k)] − (RΓ_{u})^{T} Δũ_{c}(k)}

Then

J(k) = [ṽ(k) − ṽ^{0}(k)]^{T} S [ṽ(k) − ṽ^{0}(k)] + κ  (34)

for

κ = [ỹ_{c}(k) − s̃(k)]^{T} Q [ỹ_{c}(k) − s̃(k)] + Δũ_{c}(k)^{T} RΔũ_{c}(k) − ṽ^{0}(k)^{T} Sṽ^{0}(k)

Then the cost to be minimized over ṽ(k) may be written as

J′(k) = ∥S^{\frac{1}{2}} [ṽ(k) − ṽ^{0}(k)]∥^{2}  (35)

2.4 Constraints

In all practical control problems, constraints exist on the actuator limits, the actuator rates and the system outputs. These constraints are usually posed in the form:

Δu_{min}(k + i) ≤ Δu(k + i) ≤ Δu_{max}(k + i)  i = 0, …, q − 1
u_{min}(k + i) ≤ u(k + i) ≤ u_{max}(k + i)  i = 0, …, q − 1
y_{min}(k + i) ≤ y(k + i) ≤ y_{max}(k + i)  i = 1, …, p  (36)

A more general representation of the constraints over the control horizon and the prediction horizon allowing for linear combinations is given by

LΔũ(k) ≤ l(k)
Mũ(k) ≤ m(k)
Nỹ(k) ≤ n(k)  (37)

where, like the control increments, the future control inputs are given across the control horizon as

ũ(k) = [ u(k)^{T} … u(k + q − 1)^{T} ]^{T}  (38)

and l(k) ∈ ^{λ}, m(k) ∈ ^{μ}, and n(k) ∈ ^{ν} are the constraint limits, which may vary with time. M ∈ ^{λ×qξ} and N ∈ ^{μ×pη} are real matrices and should not be confused with notation used earlier for operators. Here λ represents the total number of rate constraints across the control horizon. Similarly μ is the number of input constraints over q, and ν is the number of constraints on the outputs across the prediction horizon. In the sequel it will be desirable to “normalize” the constraints as

LΔũ(k) − l(k) + ℓ^{λ} ≤ ℓ^{λ}
Mũ(k) − m(k) + ℓ^{μ} ≤ ℓ^{μ}
N ỹ(k) − n(k) + ℓ^{ν} ≤ ℓ^{ν}  (39)

ℓ^{λ}, ℓ^{μ} and ℓ^{ν} are each a vector of ones of the same dimension as l(k), m(k) and n(k), respectively. This particular form of normalized constraint equations avoids the difficulty of normalizing in the presence of a zero constraint limit. In order to write the constraint equations in terms of the reference signal ṽ(k), it is necessary to note that

u(k + l) = ∑_{j=0}^{l} Δu(k + l) + u(k − 1)

= ∑_{j=0}^{l} ∑_{τ=0}^{n+1} T_{u}^{τ} v(k + l − τ) + u(k − 1) l = 0 … q − 1  (40)

Thus it is possible to describe the future control inputs ũ(k) in the form of equations (20) and (21):

ũ(k) = Γ_{a}ṽ(k) + Ξ_{a}ṽ^{b}(k) + X_{a}ṽ^{∞} + I_{q}u(k − 1)  (41)

where

I_{q} = [ I_{ξ} … I_{ξ} ]^{T}

Γ_{a} ∈ ^{qξ×rξ} and X_{a} ∈ ^{qξ×(p−r)ξ} are Toeplitz matrices in T_{u}^{τ}, and Ξ_{u} ∈ ^{qξ×nξ} is a Hankel matrix in T_{u}^{τ}, each of a form similar to (22) - (24). Using the constraint equations (39) in conjunction with the prediction equations for ỹ(k), Δũ(k) and ũ(k), the constraints may be combined into a single expression in terms of system signals at time k:

Dṽ(k) − Ec(k) + ℓ ≤ ℓ  (42)

Here ℓ ∈ ^{λ+μ+ν} is again a vector with l_{i} = 1, i = 1 … λ + μ + ν. In addition,

c(k) = [ ṽ^{b}(k)^{T} ṽ^{∞}^{T} u(k − 1)^{T} l(k)^{T} m(k)^{T} n(k)^{T} ]^{T}

D and E are defined as

D = [ \begin{matrix}LΓ_{a} \\ MΓ_{u} \\ NΓ_{y} \end{matrix} ]

E = [ \begin{matrix}−LΞ_{a} & −LX_{a} & −LI_{q} & I & 0 & 0 \\ −MΞ_{u} & −MX_{u} & 0 & 0 & I & 0 \\ −NΞ_{y} & −NX_{y} & 0 & 0 & 0 & I \end{matrix} ]

A sufficient condition for satisfying each row of this inequality, and thus the constraints, is given by

∥Dṽ(k) − Ec(k) + ℓ∥_{∞} ≤ 1  (43)

We will then use this requirement in conjunction with the performance index (25) to optimize subject to constraints.

2.5 Optimization Algorithm

In order to address tracking performance as well as constraints, the following weighted mixed cost function is considered based on [13]:

J_{v}(k) = ∥[ \begin{matrix}[V^{i+1}]^{\frac{1}{2}} S^{\frac{1}{2}} [ṽ(k) − ṽ^{0}(k)] \\ [W^{i+1}] & \frac{1}{2} & [Dṽ(k) − Ec(k) + ℓ] \end{matrix} ]∥^{2}  (44)

V^{i+1} = \frac{V^{i}}{∑_{k=1}^{γ} W_{k}^{i} ∥e_{k}^{i}∥}  (45)

W^{i+1} = \frac{W_{j}^{i} ∥e_{j}^{i}∥}{∑_{k=1}^{γ} W_{k}^{i} ∥e_{k}^{i}∥}  (46)

where V ∈ ^{rξ×rξ} is a diagonal weighting matrix and W_{k} denotes the kth “block” of a diagonal weighting matrix W ∈ ^{γ×γ}:

W = diag [ W_{1} W_{2} W_{3} ]

for W_{1} ∈ ^{λ×λ}, W_{2} ∈ ^{μ×μ} and W_{3} ∈ ^{ν×ν} each corresponding to the a class of constraint equations (actuator inputs, rates, or system outputs). Note that γ = λ + μ + ν. The minimization of (44) is achieved by implementing a modified form of Lawson’s algorithm ([8], [4]), in which the weighting matrices V and W are updated iteratively using

e^{i} = Dṽ^{i} − Ec(k) + ℓ  (47)
ṽ^{i} = [D^{T} W^{i}D]^{−1} D^{T} [W^{i}Ec(k) + ℓ]  (48)

where ṽ^{i} is the minimizing solution of the square of the Euclidean norm of (43). The initial values V^{0} and W^{0} may be chosen arbitrarily as the identity matrix. The choice of (45) and (46) as weightings will minimize the ∥·∥_{∞} in equation (43) because they most penalize the element of e_{k} with the greatest error. The constraints are satisfied when ∥e^{i}∥ < 1, at which point the setpoint tracking component of the mixed cost is given precedence due to V increasing relative to W. Thus convergence to a minimum of (25) can be achieved only when the constraints are satisfied.

The closed loop implementation of the algorithm is illustrated in figure (2) for this input-output case.

[figure: y; J; s; −; l, m, n; v; M^{−1}; N; v^{b}, v^{∞}; \frac{z^{−1}}{1−z^{−1}}; Δu; (Ỹ − QÑ)^{−1}; X̃ + QM̃]

Figure 2: Stable model predictive system

In this figure, J performs the on-line optimization and produces the optimal reference signal ṽ(k). Only the first element of ṽ(k) is output to the inner loop, and then the entire signal is reoptimized at time k + 1. At each time step this procedure is repeated, and thus the controller operates over a moving horizon.

2.6 Stability under constraints

To close the discussion of stable model predictive control in an input-output format and for completeness, we will state conditions for stability of the algorithm under constraints, under the assumption that stability can be achieved in the unconstrained case (surely a necessary assumption). We will rely on assuming feasibility of the optimization. Although the details of such feasibility will not be given here, they are a straightforward result of the properties of the weights in the previous section, as shown in [13]. Since the multiple-input, multiple output stability results mirror the single-input, single-output results and are proved in a similar manner to theorems 1 and 2, these proofs will be stated briefly. The complete details for the single-input, single-output case are given in [13].

Theorem 3 Assume that the output y reaches a steady state setpoint s^{∞} within p steps (see theorem 2) and that for a given prediction, control and reference horizon, a feasible set of ṽ(k) exists for which the constraints are satisfied. Then the output will asymptotically follow any setpoint change.

Proof The feasibility assumption is necessary for existence of an optimal ṽ^{∗}(k). The algorithm minimizes J^{∗}(k) subject to constraints over the specified horizons. Because of the assumptions that the output reaches the setpoints within p steps and the control increments are zero past q steps, the horizons in J^{∗}(k) may be extended to infinity, and similarly the constraints are satisfied for all time. The rest of the proof follows directly from theorem 2.

□

Theorem 4 Assume the setpoint signal s reaches a constant steady state value s^{∞} within a finite number of time steps. Then, if a feasible input signal u exists which will drive the output y to s^{∞} within N steps, the system will be stable for r > N − α and p > r + α, where α is the highest degree of any n_{ij}(z).

Proof Because of the finite impulse response nature of the inner feedback loop, u must reach its steady state value in N − α steps. Thus u may result from the relationship u = Mv if r > N − α. The remainder of the proof follows from theorem 2.

□

3 State Space Implementation

In the following we will simply illustrate how the constrained multivariable model predictive control formulation may be implemented using state space methods with state estimation. Although much of the same notation as in the previous sections will be used, note that this does not imply that the matrices are the same, but only that they correspond to the same terms in the input-output case. A model is assumed to be available in discrete state space form, with (A, B) stabilizable and (C, A) detectable:

x(k + 1) = Ax(k) + BΔu(k)
y(k) = Cx(k) + DΔu(k)

where x(k) ∈ ^{ζ} is the system state vector. The right coprime factorization of section (2) may be realized in state space form by introducing the state feedback variable v(k) as

v(k) = Z_{r}^{−1} [Δu(k) − F_{r}x(k)]

where F_{r} is chosen such that A + BF_{r} is stable and Z_{r} is any nonsingular matrix [18]. The model may then be expressed as

x(k + 1) = (A + BF_{r})x(k) + BZ_{r}v(k)  (49)
Δu(k) = F_{r}x(k) + Z_{r}v(k)  (50)
y(k) = (C + DF_{r})x(k) + DZ_{r}v(k)  (51)

From this representation, it is evident that the following relationships hold which correspond to (10) and (11):

y = Nv  (52)
Δu = Mv  (53)

where N is the transfer matrix from v to y and M is that from v to Δu. The stabilizing feedback loop of section (2) applies since the input-output properties of the loop remain the same. The controller (4) may be realized in a state space form as in [18].

3.1 State prediction

At each time k, the model of equations (49) - (51) is used to predict the plant behavior through a set of state estimates, denoted as x̂(k) ∈ ^{κ}, which will take the following form in terms of a particular point l on a prediction horizon p:

x̂(k + l + 1) = F (l)x̂(k + l) + Gv(k + l) + H(l)y(k)  l = 0, …, p − 1  (54)

Note that F (l), G(l) and H(l) are not the matrices that appear in the open-loop plant model, but are derived from these in a manner depending upon the estimation scheme being used. From earlier definitions and using (54), the vector of future state predictions may then be written as

[ x̂(k + 1) ; ⋮ ; x̂(k + p) ] = Fx̂(k) + Gṽ(k) + G^{∞}ṽ^{∞}(k) + Hy(k)  (55)

where, using the notation ∏_{j=m}^{n} F (j) to denote F (m)F (m − 1) … F (n) for m > n [9], the matrices F, G and H are given by

F = [ F (0) ; F (1)F (0) ; ⋮ ; ∏_{j=p−1}^{0} F (j) ]} p

G = \overbrace{[ \begin{matrix}G & 0 & ⋯ & 0 \\ F (1)G & G & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ [∏_{j=r}^{1} F (j)] G & [∏_{j=r}^{2} F (j)] G & ⋯ & G \\ ⋮ & ⋮ & ⋮ & ⋮ \\ [∏_{j=p−1}^{1} F (j)] G & [∏_{j=p−1}^{2} F (j)] G & ⋯ & [∏_{j=p−1}^{r+1} F (j)] G \end{matrix} ]\} p}^{r}

G^{∞} = \overbrace{[ \begin{matrix}0 & 0 & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ G & 0 & ⋯ & 0 \\ F (r + 1)G & G & ⋯ & 0 \\ [∏_{j=r+2}^{r+1} F (j)] G & F (r + 2)G & ⋯ & G \\ ⋮ & ⋮ & ⋮ & ⋮ \\ [∏_{j=p−1}^{r+1} F (j)] G & [∏_{j=p−1}^{r+2} F (j)] G & ⋯ & G \end{matrix} ]\} p}^{p−r}

H = [ H(0) ; F (1)H(0) ; ⋮ ; [∏_{j=p−1}^{1} F (j)] H(0) ]} p  (56)

To write prediction equations for the control increments Δu(k), we make use of equation (40) as well as the model equation (50) and the state predictions (54) to obtain a new expression for ũ(k):

Δũ(k) = F_{Δ}x̂(k) + G_{Δ}ṽ(k) + G_{Δ}^{∞}ṽ^{∞} + H_{Δ}y(k)  (57)

where F_{Δ}, G_{Δ}, G_{Δ}^{∞} and H_{Δ} are derived in a similar manner to equations (56):

F_{Δ} = [ F_{r} ; F_{r}F (0) ; ⋮ ; F_{r} [∏_{j=q−2}^{0} F (j)] ]} q

G_{Δ} = \overbrace{[ \begin{matrix}Z_{r} & 0 & ⋯ & 0 \\ F_{r}GZ_{r} & Z_{r} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ F_{r} [∏_{j=r−1}^{1} F (j)] GZ_{r} & F_{r} [∏_{j=r−1}^{2} F (j)] GZ_{r} & ⋯ & F_{r}GZ_{r} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ F_{r} [∏_{j=q−2}^{1} F (j)] GZ_{r} & F_{r} [∏_{j=q−2}^{2} F (j)] GZ_{r} & ⋯ & F_{r} [∏_{j=q−2}^{r} F (j)] GZ_{r} \end{matrix} ]\} q}^{r}

G_{Δ}^{∞} = \overbrace{[ \begin{matrix}0 & 0 & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ F_{r}GZ_{r} & Z_{r} & ⋯ & 0 \\ F_{r}F (r + 1)GZ_{r} & F_{r}GZ_{r} & ⋯ & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ F_{r} [∏_{j=q−2}^{r+1} F (j)] GZ_{r} & F_{r} [∏_{j=q−2}^{r+2} F (j)] GZ_{r} & ⋯ & F_{r}GZ_{r} \end{matrix} ]\} q}^{q−r}

H_{Δ} = [ 0 ; F_{r}H(0) ; ⋮ ; F_{r} [∏_{j=q−2}^{1} F (j)] H(0) ]} p  (58)

Similarly, the actual control inputs may be written across the control horizon as

ũ(k) = F_{u}x̂(k) + G_{u}ṽ(k) + G_{u}^{∞}ṽ^{∞} + H_{u}y(k) + I_{κ}u(k − 1)  (59)

where

F_{u} = [ F_{r} ; F_{r} [F (0) + I] ; ⋮ ; F_{r} [∏_{j=q−2}^{0} F (j) + ∏_{q−3}^{0} F (j)] ]} q

G_{u} = [ \begin{matrix}Z_{r} & 0 \\ [F_{r}G + I] Z_{r} & Z_{r} \\ ⋮ & ⋮ \\ F_{r} [∏_{j=r−1}^{1} F (j) + ∏_{j=r−2}^{1} F (j)] GZ_{r} & F_{r} [∏_{j=r−1}^{2} + ∏_{j=r−2}^{2} F (j)] GZ_{r} \\ ⋮ & ⋮ \\ F_{r} [∏_{j=q−2}^{1} F (j) + ∏_{j=q−3}^{1} F (j)] GZ_{r} & F_{r} [∏_{j=q−2}^{2} F (j) + ∏_{j=q−3}^{2}] GZ_{r} \end{matrix}

\begin{matrix}⋯ & 0 \\ ⋯ & 0 \\ ... & ⋮ \\ ⋯ & [F_{r}G + I] Z_{r} \\ ⋮ & ⋮ \\ ⋯ & F_{r} [∏_{j=q−2}^{r} F (j) + ∏ j = q − 3^{r}] GZ_{r} \end{matrix} ]\} q

G_{u}^{∞} = [ \begin{matrix}0 & 0 \\ ⋮ & ⋮ \\ [F_{r}G + I] Z_{r} & 0 \\ F_{r} [F (r + 1) + I] GZ_{r} & [F_{r}G + I] Z_{r} \\ ⋮ & ⋮ \\ F_{r} [∏_{j=q−2}^{r+1} F (j) + ∏_{j=q−3}^{r+1} F (j)] GZ_{r} & F_{r} [∏_{j=q−2}^{r+2} F (j) + ∏_{j=q−3}^{r+2} F (j)] GZ_{r} \end{matrix}

\begin{matrix}⋯ & 0 \\ ... & ⋮ \\ ⋯ & 0 \\ ⋯ & 0 \\ ⋮ & ⋮ \\ ⋯ & F_{r}GZ_{r} \end{matrix} ]\} q

H_{u} = [ 0 ; F_{r}H(0) ; ⋮ ; F_{r} [∏_{j=q−2}^{1} F (j) + ∏ j = q − 3^{1}F (j)] H(0) ]} p

I_{κ} = \overbrace{[ I_{ξ} I_{ξ} … I_{ξ} ]^{T}}^{q}  (60)

3.2 Least squares cost formulation

The performance index which addresses output deviation from the setpoints as well as penalizes control power is given by:

J(k) = ∑_{l=1}^{p} ∥Ĉx̂(k + l) − s(k + l)∥_{Q(l)}^{2} + ∑_{l=0}^{q} ∥Δu(k + l)∥_{R(l)}^{2}  (61)

where Ĉ is the output matrix associated with the state estimates, and Q(l) > 0 and R(l) ≥ 0 are again weighting matrices. Following the procedure of section(2.3), the performance index is written in terms of the reference signal ṽ(k). Define Q(l) and R(l) as before in (31) and let

C = diag [ Ĉ … Ĉ ]

Also define x̃_{c}(k) and Δũ_{c}(k) as the contributions to [ x̂(k + 1) … x̂(k + p) ]^{T} and Δũ(k) that are independent of the future reference signals ṽ(k):

x̃_{c}(k) = Fx̂(k)G^{∞}ṽ^{∞}(k) + Hy(k)
Δũ_{c}(k) = F_{Δ}x̂(k) + G_{Δ}^{∞}ṽ^{∞} + H_{Δ}y(k) + I_{κ}u(k − 1)

Then

J(k) = ṽ(k)^{T} Sṽ(k) + 2 {[x̃_{c}(k) − s̃(k)]^{T} QCG + Δũ_{c}(k)^{T} RG_{·}} ṽ(k) +
[x̃_{c}(k) − s̃(k)]^{T} C^{T} QC [x̃_{c}(k) − s̃(k)] + Δũ_{c}(k)^{T} RΔũ_{c}(k)  (62)

where

S = G^{T} C^{T} QCG + G_{Δ}^{T}RG_{Δ}

To complete the squares, let

ṽ^{0}(k) = S^{−1} {(QCG)^{T} [s̃(k) − x̃_{c}(k)] − (RG_{Δ})^{T} Δũ_{c}(k)}

Then again

J(k) = [ṽ(k) − ṽ^{0}(k)]^{T} S [ṽ(k) − ṽ^{0}(k)] + κ  (63)

for

κ = [x̃_{c}(k) − s̃(k)]^{T} C^{T} QC [ỹ_{c}(k) − s̃(k)] + Δũ_{c}(k)^{T} RΔũ_{c}(k) − ṽ^{0}(k)^{T} Sṽ^{0}(k)

Then the cost to be minimized over ṽ(k) may again be written in terms of a Euclidean norm as

J′(k) = ∥S^{\frac{1}{2}} [ṽ(k) − ṽ^{0}(k)]∥^{2}  (64)

3.3 Constraints

In parallel with section (2.4), constraints and linear combinations of constraints may be expressed as

LΔũ(k) ≤ l(k)
Mũ(k) ≤ m(k)
N [ x̂(k + 1) ; ⋮ ; x̂(k + p) ] ≤ n(k)  (65)

where l(k) ∈ ^{λ}, m(k) ∈ ^{μ}, and n(k) ∈ ^{ν} may vary with time. Again λ represents the total number of rate constraints across the control horizon, μ is the number of input constraints over q, and ν is the number of constraints on the state estimates across the prediction horizon. Again “normalizing” the constraints:

LΔũ(k) − l(k) + ℓ^{λ} ≤ ℓ^{λ}
Mũ(k) − m(k) + ℓ^{μ} ≤ ℓ^{μ}
N [ x̂(k + 1) ; ⋮ ; x̂(k + p) ] − n(k) + ℓ^{ν} ≤ ℓ^{ν}  (66)

ℓ^{λ}, ℓ^{μ} and ℓ^{ν} are each a vector of ones of the same dimension as l(k), m(k) and n(k), respectively. This form of normalized constraint equations again avoids the difficulty of normalizing on a zero constraint limit.

Combining the constraint equations (66) with the prediction equations (55), (57) and (59), the constraints may be written in terms of current signals at time (k) as

Dṽ(k) − Ec(k) + ℓ ≤ ℓ  (67)

where

c(k) = [ x̂(k)^{T} y(k)^{T} v^{∞}(k)^{T} u(k − 1)^{T} l(k)^{T} m(k)^{T} n(k)^{T} ]^{T}

and ℓ ∈ ^{λ+μ+ν} is again a vector of ones. D and E are defined as

D = [ \begin{matrix}LG_{Δ} \\ MG_{u} \\ NG \end{matrix} ]  (68)

E = [ \begin{matrix}−LF_{Δ} & −LH_{Δ} & −LG_{Δ}^{∞} & −LI_{κ} & I & 0 & 0 \\ −MF_{u} & −MH_{u} & −MG_{u}^{∞} & 0 & 0 & I & 0 \\ −NF & −NH & −NG^{∞} & 0 & 0 & 0 & I \end{matrix} ]  (69)

A sufficient condition for satisfying the constraint equations is then given by

∥Dṽ(k) − Ec(k) + ℓ∥_{∞} ≤ 1  (70)

3.4 Optimization and implemention

The mixed cost function for the state space case is of an identical form to that of section (2.5), and, noting the corresponding notation, is implemented in exactly the same iterative fashion using the modified Lawson algorithm, so it will not be repeated here. Instead, we will illustrate the practical implementation of the closed loop system, which is shown in figure 2:

[figure: y; P; J; s; −; x̂; l, m, n; v; M^{−1}; N; v^{∞}; v^{∞}; \frac{z^{−1}}{1−z^{−1}}; Δu; (Ỹ − QÑ)^{−1}; X̃ + QM̃]

Figure 3: Stable model predictive system

Here J performs the optimization algorithm and outputs the optimal reference signal v(k). Note that this implies that only the first increment of ṽ(k) is actually used at time k, and the entire signal is reoptimized at time k + 1, as in the input-output case. P performs the state prediction. The procedure is then repeated at each time step over a moving horizon.

4 Conclusions

In this paper we have addressed the relatively new method of constrained generalized predictive control, and extended it to include the multivariable case with a state space implementation. As a result, it is possible to stabilize multivariable systems using on-line optimization wrapped around a stabilizing inner feedback loop. Furthermore, this stability is guaranteed for given prediction, control and reference horizons through a monotonically decreasing cost property. Although the stability theorems are given in an input-output format, it is anticipated that they extend to the state space framework.

References

[1] D.W. Clarke and C. Mohtadi. Properties of generalized predictive control. Automatica, 25(6):859–875, 1989.

[2] D.W. Clarke, C. Mohtadi, and P.S. Tuffs. Generalized predictive control - part I. The basic algorithm. Automatica, 23(2):137–148, 1987.

[3] D.W. Clarke, C. Mohtadi, and P.S. Tuffs. Generalized predictive control - part II. Extensions and interpretations. Automatica, 23(2):149–160, 1987.

[4] R .R. Fletcher. Practical Optimization. Wiley, first edition, 1980.

[5] B.A. Francis. A Course in H_{∞} Control Theory. Springer-Verlag, first edition, 1987.

[6] S.A. Heise and J.M. Maciejowski. Stability of constrained MBPC using an internal model control structure. In Advances in Model-Based Predictive Control. Oxford University Press, 1993.

[7] B. Kouvaritakis, J.A. Rossiter, and A.O.T. Chang. Stable generalized predictive control: an algorithm with guaranteed stability. IEE Proceedings-D, 139(4):349–362, 1993.

[8] C.L. Lawson. Contributions to the theory of linear least maximum approximations. Master’s thesis, University of California, Los Angeles, 1961.

[9] J.M. Maciejowski and S.A. Heise. Heuristic robustness analysis of model-based predictive controllers. In Proceedings of the 12th World Congress. International Federation of Automatic Control, 1993.

[10] M. Morari. Model predictive control: Multivariable control technique of choice in the 1990’s? 1993.

[11] K.R. Muske and J.B. Rawlings. Model predictive control with linear models. Process Systems Engineering, 39(2):262–287, 1993.

[12] J.A. Rossiter and B. Kouvaritakis. Advances in generalized and constrained generalized predictive control with application to a dynamometer model. In Advances in model based predictive control. Oxford University Press, 1993.

[13] J.A. Rossiter and B. Kouvaritakis. Constrained stable generalized predictive control. IEE Proceedings-D, 140(4):243–254, 1993.

[14] T.T.C. Tsang and D.W. Clarke. Generalized predictive control with input constraints. IEE Proceedings, 135, Pt. D(6):1451–460, 1988.

[15] D.C. Youla, J.J. Bongiorno, and H.A. Jabr. Modern Wiener-Hopf design of optimal controllers - part II. the multivariable case. IEEE Transactions on Automatic Control, 21:319–338, 1976.

[16] E. Zafiriou. Robust model predictive control of processes with hard constraints. Computers Chem. Engineering, 14(4/5):359–371, 1990.

[17] A. Zheng and M. Morari. Stability of model predictive control with soft constraints. Submitted to the IEEE Transactions on Automatic Control, 1993.

[18] K. Zhou, J.C. Doyle, and K. Glover. Robust and Optimal Control. unpublished draft, 1993.