Heuristic Robustness Analysis of Model-Based Predictive Controllers

J.M.Maciejowski ^{∗}  S.A.Heise

Cambridge University Engineering Department Cambridge CB2 1PZ, England. Fax: +44 223 332662. Telex: 81239 DEPENG G Email: jmm@eng.cam.ac.uk, Tel: +44 223 332732

January 25, 1993

Submitted to the 12th IFAC World Congress, Sydney, 1993, as a Regular

Paper, subject area R-15. Preferred mode of presentation: 20 minute oral.

Keywords Predictive control, Constraints, Robustness.

Abstract

Model-based predictive control schemes have been applied very successfully in the process industries, even though very little is known about them theoretically — typically even stability cannot be proved. Their use in reconfigurable control schemes makes them of potential interest in applications such as flight control, and the increasing speed of hardware makes such applications feasible in the near future. However an essential prerequisite for such applications is a much better theoretical understanding of MBPC systems. The main theoretical contribution of this paper is a proof that most MBPC schemes are piecewise-LTI, even when constraints are active. On the basis of this, some heuristics are indicated for the selection of competing MBPC designs.

∗This author for correspondence.

1 Introduction

Model-based predictive controllers (MBPC) have proved very successful in a variety of process control applications [6, 12]. These controllers typically have a number of parameters, including integer-valued parameters such as prediction horizons, which are selected according to heuristic guidelines [5, 15]. The stability of such controllers has not been established, except under assumptions which are either very conservative or difficult to verify [19, 13].

One of the attractions of model-based predictive control is that it offers a systematic way of approaching the problem of reconfiguring the control law in the face of equipment failure, or of a severe change in the operating environment. This makes it an intriguing possibility for applications such as flight control or automotive control.

The keys to this capability are, firstly, that an explicit internal model of the controlled plant is used by the controller, and secondly, that the control signal is obtained by solving a constrained optimization problem. The increasing use of accurate models, in both process and flight control [9, 11], together with the emergence of good fault detection algorithms, makes it possible to envisage fast and accurate identification of most failures or damage which might occur, with consequent update of an internal model to represent the condition of the plant. The constrained optmization gives a systematic method of finding control signals which drive the plant away from dangerous regions of the state-space, or shut it down safely, when a failure or damage does occur.

The increasing speed of computers makes it possible to solve even quite large constrained optmization problems in real time at ever higher sample rates, so that it is quite possible to foresee this approach being feasible in high-bandwidth applications such as flight and automotive control (at sample rates of 10 Hz, say). The use of application-specific integrated circuits will probably make even higher rates feasible.

In such applications, however, conventional control schemes are generally quite adequate in normal circumstances, so the only advantage to be gained from MBPC schemes is their performance in severely adverse conditions. In process applications, by contrast, predictive control is used largely for its economic benefits in normal circumstances, and there is usually enough time to switch to more conventional control if things go badly wrong. It is therefore more essential to establish the robustness of a predictive control scheme in one of these new applications than it would be in process control (although even there it is clearly desirable to do so).

In this paper we propose that the analysis of predictive control schemes can be assisted by the use of well-established system-theoretic computations. To do this we exploit the fact that most of the proposed predictive control schemes operate as linear time-invariant control laws most of the time. We do not have any real theorems yet, but instead suggest guidelines for design which supplement existing guidelines. Some similar work has been reported in [14].

2 Concepts of robustness for predictive control

Relatively little work has appeared on the stability of constrained MBPC control laws, and even less on their robustness. Prett and Garcia [15] correctly state that even an appropriate notion of stability is difficult to define for MBPC systems. Robustness (which in some sense includes stability) is even more elusive.

It is frequently stated that MBPC laws which are stable in the absence of constraints remain stable if control signal amplitude constraints are imposed (eg [15]). The argument usually emplyed is that input constraints can only reduce the energy in the control signals, and can therefore not excite instability in a previously stable system. In these discussions stability is not defined, and the following examples show that the claims are probably not true. First, consider an integrator with proportional output feedback and an input disturbance:

ẋ(t) = ax(t) + d(t)  (1)

This is asymptotically stable providing that a < 0. However, if the control signal is limited to amplitude no greater than 1, so that the equation becomes

ẋ = sat(ax(t)) + d(t)  (2)

and if the input disturbance has a constant value greater than 1:

d(t) = d_{o} (|d_{o}| > 1)  (3)

then |x(t)| → ∞ as t → ∞. The constrained system is therefore not BIBO stable, although the equilibrium remains globally asymptotically stable, in the sense of Lyapunov, if d(t) ≡ 0. Secondly, consider a SISO LTI feedback system which is conditionally stable — namely, one whose Nyquist locus crosses the negative real axis to the left of the point −1. Imposing an input constraint on such a system makes the presence of stable limit cycles possible. Is the constrained system to be regarded as stable or not? Certainly global stability of the equilibrium point is lost when the constraint is present.

These examples simply make the rather obvious point that once constraints are applied, even input constraints, then several nonequivalent notions of stability are possible, and one has to be careful about choosing appropriate ones.

A (rather imprecise) definition of stability which may be appropriate for MBPC is the following, stated by Kohn [10]: ‘A system is stable if, for any time period of operation and for all possible environmental conditions, the absorbed energy is bounded by a finite quantity, and bounded commanded actions lead to bounded output data’. Note that according to this definition the integrator with constrained input is stable, since it distingushes between ‘commanded’ actions and those which might result from disturbance signals.

Precise and conventional notions of stability are used in the work on the closely-related subject of receding-horizon control, in particular in [13] and [18]. Both of these use the conventional notion of Lyapunov stability of a unique equilibrium point, and assume that the constraints and environment are such that this can be attained.

3 MBPC control is ‘piecewise-LTI’

We consider the following version of MBPC. The plant is modelled as a set of linear time-invariant state-space systems, each system being accomanied by one or more (possibly empty) set of linear inequality constraints on the inputs, the rates of change of the inputs, and on the states. Each member of the set represents the plant near a particular operating point, and in a particular condition, and switching between elements of this set is assumed to take place relatively rarely. Note that the combination of linear system with linear inequalities allows the modelling of some very important classes of nonlinear systems, such as those which are linear except for limits on actuator rates and extensions.

The objectives of control are represented by a set of quadratic functionals of the predicted input and state trajectories over a (usually finite) prediction horizon, together with further (possibly empty) sets of linear inequalities on the states. Here again each combination of a quadratic functional and of a set of constraints is in effect for a given operating region and plant condition, so that switching between these is again assumed to be relatively infrequent.

The control signal is determined by minimising the particular quadratic functional which is in effect at that time, subject to all those constraints which are then in force. This minimization gives a trajectory of future control signals over some prespecified optimization horizon; some initial portion of this trajectory is applied to the plant, and the procedure is repeated to determine the next portion of the control trajectory.

This formulation includes the DMC [15] formulation, its extensions such as [17], the MPHC formulation [16], and GPC [4], although it is more general than any of these. If linear equality constraints are also allowed then it includes those ‘receding-horizon’ control laws in which a terminal constraint on the state is included (eg [13]), when they are applied to piecewise-LTI plants.

In this formulation the control law is linear time-invariant, so long as the set of active constraints remains unchanged. The control law is therefore a ‘variable structure’ law, which switches between a finite number of LTI laws. It is well known that the MBPC law is LTI when no constraints are active [3, 15, 1]; the fact that it remains piecewise-LTI when constraints are active is apparently known to some workers in the area, but does not appear in the literature, to our knowledge. We shall now prove it for the case when the problem is formulated in discrete time, and only the first of the computed control signals is implemented before recalculating future control signals.

3.1 State space modelling and state prediction

We assume a set of system models is available in discrete state space form:

x_{i}(k + 1) = A_{i}x_{i}(k) + B_{u_{i}}u_{i}(k) + B_{d_{i}}d_{i}(k)  (4)
y_{i}(k) = C_{i}x_{i}(k) + w_{i}(k − 1) + v_{i}(k)  (5)

where i is used to index the plant model in use, x_{i}(k) ∈ ℜ^{κ} is the system state vector, u_{i}(k) ∈ ℜ^{ξ} is the system manipulated variable (control input) set, and y_{i}(k) ∈ ℜ^{η} is a vector of controlled outputs, measured at time k. Note that the input and measurement vectors are indexed by i as well, to indicate that the variables which appear in these vectors may depend upon which model is in use. d_{i}(k) is an unmeasured state disturbance input and w_{i}(k − 1) is an unmeasured output disturbance vector, possibly representing error due to plantmodel mismatch or other unmodelled effects. v_{i}(k) is an additional zero mean white noise vector corrupting the measured outputs.

In order to cast the state space model-based control problem into a conventional quadratic programming form, we seek to express the state vector in terms of the change in the control vector at time k, Δu_{i}(k) = u_{i}(k) − u_{i}(k − 1). A convenient form is given by

[ \begin{matrix}Δx_{i}(k + 1) \\ y_{x_{i}}(k + 1) \end{matrix} ] = [ \begin{matrix}A_{i} & 0 \\ C_{i}A_{i} & I \end{matrix} ] [ \begin{matrix}Δx_{i}(k) \\ y_{x_{i}}(k) \end{matrix} ] + [ \begin{matrix}B_{u_{i}} \\ C_{i}B_{u_{i}} \end{matrix} ] Δu_{i}(k)
+ [ \begin{matrix}B_{d_{i}} & 0 \\ C_{i}B_{d_{i}} & I \end{matrix} ] [ \begin{matrix}Δd_{i}(k) \\ Δw_{i}(k) \end{matrix} ] (6)

where y_{x_{i}} corresponds to a vector of outputs free from measurement noise. Certain elements of the state vector may not be available for measurement at time k, so we will wish to determine an estimate of the state using some sort of observer. To this end, the disturbance inputs d_{i}(k) and w_{i}(k) may be modelled in either a deterministic or a stochastic sense, depending upon the particular application under consideration. In MBPC these disturbances are normally assumed to be in the form of impulses, steps, ramps or the output of some other deterministic model suitable for representing slowly varying load disturbances on a process. Since it is more natural to pose the estimation problem in a stochastic setting, however, we model the disturbances as linear systems driven by white noises Δd̃_{i}(k) and Δw̃_{i}(k) and augment these to the state equations as follows:

[ \begin{matrix}Δx_{i}(k + 1) \\ y_{x_{i}}(k + 1) \\ Δd_{i}(k + 1) \\ Δw_{i}(k + 1) \end{matrix} ] = [ \begin{matrix}A_{i} & 0 & B_{d_{i}} & 0 \\ C_{i}A_{i} & I & C_{i}B_{d_{i}} & I \\ 0 & 0 & A_{d_{i}} & 0 \\ 0 & 0 & 0 & A_{w_{i}} \end{matrix} ] [ \begin{matrix}Δx_{i}(k) \\ y_{x_{i}}(k) \\ Δd_{i}(k) \\ Δw_{i}(k) \end{matrix} ]

+ [ \begin{matrix}B_{u_{i}} \\ C_{i}B_{u_{i}} \\ 0 \\ 0 \end{matrix} ] Δu_{i}(k) + [ \begin{matrix}0 & 0 \\ 0 & 0 \\ I & 0 \\ 0 & I \end{matrix} ] [ \begin{matrix}Δd̃_{i}(k) \\ Δw̃_{i}(k) \end{matrix} ] (7)

y_{i}(k) = [ 0 I 0 0 ] [ \begin{matrix}Δx_{i}(k) \\ y_{x_{i}}(k) \\ Δd_{i}(k) \\ Δw_{i}(k) \end{matrix} ] + v_{i}(k)  (8)

where the corresponding covariance matrices are given as

E \{[ \begin{matrix}Δd̃_{i} \\ Δw̃_{i} \end{matrix} ] [ Δd̃_{i}^{T} Δw̃_{i}^{T} ]\} = [ \begin{matrix}Q_{d} & 0 \\ 0 & Q_{w} \end{matrix} ]  (9)
and E {v_{i}(k)v_{i}(k)^{T} } = Q_{v}  (10)

Here Q_{d}, Q_{w} and Q_{v} will be tuning parameters of the observer.

If the system being considered is such that the disturbances d_{i}(k) and w_{i}(k) are stochastic inputs such as white noise or zero mean, random amplitude impulses, steps or ramps, then it is not necessary to augment them to the state vector [ Δx_{i}(k)^{T} y_{x_{i}}(k)^{T} ]^{T}, and the obvious simplifications occur.

At time k each member of the model set used to represent the plant is used to predict the plant behaviour in the following way:

x̂_{i}(k + l + 1) = F_{i}(l)x̂_{i}(k + l) + G_{i}Δu_{i}(k + l) + H_{i}(l)y_{i}(k) (l = 0, …, p − 1) (11)

Here x̂_{i}(k) is an estimate of the state vector of the plant at time k, obtained using the i’th model. p is the prediction horizon. We should emphasise that F_{i}(l), G_{i} and H_{i}(l) are not the matrices which appear in the open-loop plant model itself, but are derived from these in a manner which depends on the estimation/prediction scheme being used. This depends in general on the kind of observer chosen, on the assumptions about disturbances which affect the plant, and on the noise which affects the measurements. We need to assume that F_{i}(l), G_{i} and H_{i}(l) depend only on i and l, but not on k. This excludes certain possibilities, such as non-stationary Kalman filters, but includes all the schemes which are commonly employed with MBPC. Normally

H_{i}(l) = 0 for l > 0  (12)

and we shall assume that this is the case. It is also commonly the case that F_{i}(0) ≠ F_{i}(1) and F_{i}(l) = F_{i}(l + 1) for l > 0, but we shall not assume this. The state respresentation of the plant is not assumed to be minimal, so that DMC can be included by letting x_{i}(k) be a vector of past input changes.

In the following we shall use the notation ∏_{j=m}^{n} F_{i}(j) to denote F_{i}(m)F_{i}(m− 1) … F_{i}(n), where m > n. We now introduce the matrices

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

G_{i} = [ \begin{matrix}G_{i} & 0 & ⋯ & 0 \\ F_{i}(1)G_{i} & G_{i} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ [∏_{j=q}^{1} F_{i}(j)] G_{i} & [∏_{j=q}^{2} F_{i}(j)] G_{i} & ⋯ & G_{i} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ [∏_{j=p−1}^{1} F_{i}(j)] G_{i} & [∏_{j=p−1}^{2} F_{i}(j)] G_{i} & ⋯ & [∏_{j=p−1}^{q+1} F_{i}(j)] G_{i} \end{matrix} ] (14)

H_{i} = [ H_{i}(0) ; F_{i}(1)H_{i}(0) ; ⋮ ; [∏_{j=p−1}^{1} F_{i}(j)] H_{i}(0) ]  (15)

where q ≤ p. q is going to be the optimization horizon, and in the ensuing optimization it will be assumed that

Δu_{i}(k + l) = 0 for l > q.  (16)

The vector of q + 1 future control changes will be denoted by

ΔU_{i}(k) = [Δu_{i}(k)^{T} , …, Δu_{i}(k + q)^{T}]^{T} .  (17)

With this notation and assumption we can write the vector of future state predictions as

[ x̂_{i}(k + 1) ; ⋮ ; x̂_{i}(k + p) ] = F_{i}x̂_{i}(k) + G_{i}ΔU_{i}(k) + H_{i}y_{i}(k).  (18)

3.2 Constraints

The i’th constraint set is assumed to consist of linear constraints on the changes in the control signals (corresponding to actuator rate constraints, typically):

L_{i}ΔU_{i}(k) ≥ [ l_{i1}(k) ; ⋮ ; l_{iλ}(k) ]  (19)

of linear constraints on the control signals themselves (corresponding to actuator authority limits, typically):

M_{i} [ u_{i}(k) ; ⋮ ; u_{i}(k + q) ] ≥ [ m_{i1}(k) ; ⋮ ; m_{iμ}(k) ]  (20)

and of linear state constraints (which includes output constraints, of course):

N_{i} [ x̂_{i}(k + 1) ; ⋮ ; x̂_{i}(k + p) ] ≥ [ n_{i1}(k) ; ⋮ ; n_{iν}(k) ]  (21)

where the right hand side of the last inequality may depend (linearly) on predicted set-points for state and/or output variables. Note that in these inequalities λ ≠ μ ≠ ν in general.

Now these three sets of inequalities can be rewritten as a single linear inequality constraint on the vector of future input changes ΔU_{i}(k), with the right hand side being linear in the measurements of inputs and outputs. (Note that the control u(k) is still to be determined at time k.) To show this for the control signal level constraints we just need to note that

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

so that the control level constraints become

M_{i} [ \begin{matrix}I & 0 & ⋯ & 0 \\ I & I & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ I & I & ⋯ & I \end{matrix} ] ΔU_{i}(k) ≥ [ m_{i1}(k) ; m_{i2}(k) ; ⋮ ; m_{iμ}(k) ] − M_{i} [ u_{i}(k − 1) ; u_{i}(k − 1) ; ⋮ ; u_{i}(k − 1) ] .(23)

The state constraints can clearly be expressed as

N_{i}G_{i}ΔU_{i}(k) ≥ [ n_{i1}(k) ; n_{i2}(k) ; ⋮ ; n_{iν}(k) ] − N_{i}F_{i}x̂_{i}(k) − N_{i}H_{i}y_{i}(k)  (24)

Now all of the inequality constraints can be written as the single constraint

[ \begin{matrix}L_{i} \\ W_{i} \\ V_{i} \end{matrix} ] ΔU_{i}(k) ≥ [ l_{i1}(k) ; ⋮ ; l_{iλ}(k) ; m_{i1}(k) ; ⋮ ; m_{iμ}(k) ; n_{i1}(k) ; ⋮ ; n_{iν}(k) ] + [ \begin{matrix}0 \\ X_{i}u_{i}(k − 1) \\ Y_{i}x̂_{i}(k) + Z_{i}y_{i}(k) \end{matrix} ]  (25)

where the matrices W_{i}, V_{i}, X_{i}, Y_{i}, Z_{i} are obtained from the previous expressions in an obvious manner, and are independent of any input or output measurements. x̂_{i}(k) itself depends on input and output measurements in a linear, time-invariant manner (according to the model introduced above), so it is seen that all the constraints can be expressed as linear inequality constraints on future control signal changes, and that the constraint levels depend in a linear, time-invariant manner on input and output measurements. Note that each row of this vector inequality corresponds to a physical constraint on the plant or its actuators.

Arguing in a similar manner, it can be seen that adding a terminal state constraint of the form

x̂_{i}(k + p) = const.  (26)

as is done in the ‘receding-horizon’ literature, results in adding a linear equality constraint on future input changes.

3.3 The cost functional

The control is determined by minimising the following cost functional:

J_{i}(k) = ∑_{l=1}^{p} ||C_{i}x̂_{i}(k + l) − s_{i}(k + l)||_{Q_{i}(l)}^{2} + ∑_{l=0}^{q} ||Δu_{i}(k + l)||_{R_{i}(l)}^{2} (27)

where s_{i}(k) is the vector of set-points for states and/or outputs at time k. Note that the weights which appear in J_{i}(k) are assumed to be independent of k, although they certainly depend on i and may depend on l. This cost function can be written as

J_{i}(k) = ΔU_{i}(k)^{T} [G_{i}^{T} C_{i}^{T} Q_{i}C_{i}G_{i} + R_{i}]ΔU_{i}(k) +
+2{[F_{i}x̂_{i}(k) + H_{i}y_{i}(k)]^{T} C_{i}^{T} − S_{i}(k)^{T} }Q_{i}C_{i}G_{i}ΔU_{i}(k) +
+K_{i}(k)  (28)

where Q_{i}, R_{i}, C_{i} are the block-diagonal matrices

Q_{i} = diag[Q_{i}(1), …, Q_{i}(p)]  (29)
R_{i} = diag[R_{i}(1), …, R_{i}(q)]  (30)
C_{i} = diag[C_{i}, …, C_{i}]  (31)

and

S_{i}(k) = [s_{i}(k + 1)^{T} , …, s_{i}(k + p)^{T}]^{T} .  (32)

The final term K_{i}(k) is independent of ΔU_{i}(k) and can therefore be omitted from J_{i}(k) without affecting the optimization. (This term is in fact given by

K_{i}(k) = S_{i}(k)^{T} Q_{i}{S_{i}(k) − 2C_{i}[F_{i}x̂_{i}(k) + H_{i}y_{i}(k)]} +
+[F_{i}x̂_{i}(k) + H_{i}y_{i}(k)]^{T} C_{i}^{T} Q_{i}C_{i}[F_{i}x̂_{i}(k) + H_{i}y_{i}(k)]  (33)

.)

At this point we need to assume that S_{i}(k) is in fact independent of k. In practice this can only mean that all the set-points are constant (at least while the index i is active). However, this is not really restrictive, because if the setpoint trajectories can be modelled as the step responses of some invertible LTI system then the inverse dynamics of that system can be included in the plant model, and constant set-points provided for the outputs of the inverse system.

It is now apparent that the optimization problem to be solved is a standard quadratic program:

min_{ΔU_{i}(k)}  ΔU_{i}(k)^{T} A_{i}ΔU_{i}(k) + [x̂_{i}(k)^{T} , y_{i}(k)^{T} , c_{i}^{T}]B_{i}ΔU_{i}(k) (34)

subject to  D_{i}ΔU_{i}(k) ≥ E_{i} [ \begin{matrix}u_{i}(k − 1) \\ x̂_{i}(k) \\ y_{i}(k) \\ c_{i} \end{matrix} ]  (35)

where the matrices A_{i}, B_{i}, D_{i} and E_{i}, and the vector c_{i}, are obtained in a straightforward manner from the various matrices introduced above. With the above assumption on the set-points, and with the assumption that the constraint levels l_{ij}(k), m_{ij}(k), n_{ij}(k) are independent of k, all of these are independent of k.

Now suppose that the controlled plant is sitting on a particular set of constraints. In particular, suppose that

D_{ij}ΔU_{i}(k) = E_{ij} [ \begin{matrix}u_{i}(k − 1) \\ x̂_{(}k) \\ y_{i}(k) \\ c_{i} \end{matrix} ]  (36)

where D_{ij} and E_{ij} contain those rows of D_{i} and E_{i} which correspond to the active constraints. Then it is well known [7] that the optimal solution ΔU_{i}^{∗}(k) is given by

[ \begin{matrix}2A_{i} & −D_{ij}^{T} \\ −D_{ij} & 0 \end{matrix} ] [ \begin{matrix}ΔU_{i}^{∗}(k) \\ λ^{∗} \end{matrix} ] = − \left[{0  B_{i}^{T} \atop E_{ij}}\right] [ \begin{matrix}u_{i}(k − 1) \\ x̂_{i}(k) \\ y_{i}(k) \\ c_{i} \end{matrix} ] (37)

where λ^{∗} is the vector of Lagrange multipliers at the optimal solution. Now, defining

[ \begin{matrix}2A_{i} & −D_{ij}^{T} \\ −D_{ij} & 0 \end{matrix} ]^{−1} = [ \begin{matrix}L_{i1} \\ L_{i2} \end{matrix} ]  (38)

(the inverse exists if the QP problem has a unique feasible solution) we obtain the future control vector as

ΔU_{i}^{∗}(k) = −L_{i1} \left[{0  B_{i}^{T} \atop E_{ij}}\right] [ \begin{matrix}u_{i}(k − 1) \\ x̂_{i}(k) \\ y_{i}(k) \\ c_{i} \end{matrix} ]  (39)

If the control changes Δu_{i}^{∗}(k), Δu_{i}^{∗}(k + 1), …, Δu_{i}^{∗}(k + r − 1) are applied, and the optimization repeated at time k+r then, under the stated assumptions, this is a linear periodic feedback law, with period r. In the most usual case, which is r = 1, the feedback law is linear and time-invariant. This shows, as we claimed above, that the MBPC laws usually considered are ‘piecewise linear time-invariant’, with switches between particular LTI feedback laws occuring when the plant model changes, or the active constraints change, or the setpoints change, or the cost function changes.

4 LTI analysis of MBPC

We suggest that analysis of MBPC systems should take advantage of their piecewise-LTI nature. Precise theorems are hard to come by. For example, it is possible, at least in principle, to check the stability of an MBPC system in each possible constrained condition, using standard LTI techniques. It is clear that this does not give a sufficient condition for the stability of an MBPC system. But it is not clear whether it gives a necessary condition for stability. Nevertheless, it is plausible to believe that an MBPC system will be stable if each of its constituent LTI laws gives closed-loop stability and if switching between these laws can be guaranteed to occur rarely enough. Therefore one heuristic for obtaining stability is to choose tuning parameters such that constraints are unlikely to be activated.

This can be assessed by computing an appropriate norm of the closed-loop operator from the disturbance vector to the vector of variables which are constrained in some way. In the first place this can be done assuming that all constraints are inactive. If the closed loop is stable under this assumption (surely a necessary condition for an acceptable design) then any such norm can often be evaluated from the frequency response of the corresponding LTI operator. From the development given above it is straightforward to compute the frequency response of the controller (ie with the sequence {y_{i}(k)} as input and the sequence {u_{i}(k)} as output. From the model of the plant one can compute the frequency responses from the control inputs to the outputs, and from the disturbance inputs to the outputs. From these frequency responses one can compute any required closed-loop frequency response, including those from disturbance inputs to constraint variables. If a state-space model of the plant is available then most of the appropriate norms can be evaluated directly from the state-space models without going through frequency-response computations first [2].

The ‘appropriate’ norm depends on how much is known about typical and worst-case disturbances; for ‘text-book’ cases these are well known. If almost nothing is known about the disturbances then the ‘peak-gain’ norm may be appropriate; this is the largest l_{1} norm of any row of the (closed-loop) impulse response. The choice of a small ‘peak-gain’ norm is an attempt to avoid activating constraints, whatever form the disturbance may take. If the disturbance is persistent and can be regarded as a Gaussian stochastic process with unknown spectrum, then limiting the H_{∞} norm will limit the RMS value of the constraint variables. If the covariance of the disturbances can be estimated, then the RMS value of the constraint variables can also be estimated, and from this the probability of activating constraints can be estimated. If in addition the spectrum of the disturbances is known, then a model of the disturbances should appear as part of the plant model, so that this is driven by white disturbances. In this case the H_{2} norm gives the RMS value of the constraint variables.

Frequently much more is known about the predominant or worst-case disturbances, from a physical understanding of the plant and its environment. In particular, if disturbances are occasional deterministic signals of a known ‘shape’ then they can again be represented inside the plant model, and their occurrence triggered by an impulsive disturbance input. In this case the l_{∞} norm of the impulse response gives the peak value of any constraint variable, and is therefore an appropriate one to use.

When analysing an assumed constrained condition the desiderata are less clear. It is reasonable to try to minimise the risk of activating even more constraints than are already active, or of exchanging currently active constraints for constraints which are judged to be less desirable, for example because they represent more dangerous conditions. But in general the possibility of de-activating constraints should be encouraged. One approach to doing this is the following. Suppose, for the sake of simplicity, that some of the predicted state constraints are active;

N_{ij} [ x̂_{i}(k + 1) ; ⋮ ; x̂_{i}(k + p) ] = [ n_{iν_{1}} ; ⋮ ; n_{iν_{l}} ]  (40)

Let the vector of predicted states, evaluated at time j be X̂_{1}^{p}(j), and let ΔX̂_{1}^{p}(j, k) = X̂_{1}^{p}(j) − X̂_{1}^{p}(k). Then the current constraints will become inactive at the first j for which

N_{ij}ΔX̂_{1}^{p}(j, k) > 0.  (41)

A measure of how far away from the constraints the system will then be is given by ||(I − P)ΔX̂_{1}^{p}(j, k)||, where P is the orthogonal projector onto the kernel of N_{ij}, namely

I − P = V_{l}V_{l}^{T}  (42)

where V_{l} is the matrix of the first l right singular vectors of N_{ij} and l = rank(N_{ij}) [8]. So if G is the (closed-loop) transfer function matrix from the disturbance inputs to {ΔX̂_{1}^{p}(j, k)} (note that this transfer function is welldefined, since the corresponding operator is independent of k according to our assumptions), then one should select an MBPC scheme which gives a relatively large minimum gain of (I − P)G. It is also necessary to ensure that only N_{ij}ΔX̂_{1}^{p}(j, k) > 0 can occur, and so this approach can only be used when the disturbance is characterised sufficiently narrowly for this to be possible.

5 Conclusions

In this paper we have pointed out that most of the MBPC schemes which have been proposed and applied are piecewise-LTI even when constraints are active. We have proved this in a set up which covers schemes such as DMC, MAC and their usual generalisations, and some of the ‘receding-horizon’ schemes, and in sufficient detail to allow consequent computations to be readily deduced.

Proceeding from this observation we have suggested some heuristic guidelines for selecting competing MBPC designs. These guidelines do not replace those commonly proposed (eg [5, 15]), but rather supplement them.

The major point of our suggestions is that the piecewise-LTI nature of MBPC schemes can probably be exploited in the analysis of their stability, robustness and performance.

References

[1] Bitmead,R.R, Gevers,M, and Wertz,V, 1990, Adaptive Optimal Control: The Thinking Man’s GPC, (New York: Prentice-Hall).

[2] Boyd,S.P, and Barratt,C.H, 1991, Linear Controller Design — Limits of Performance, (Englewood Cliffs,NJ: Prentice-Hall).

[3] Chang and Seborg,D, 1983, A linear programming approach for multivariable feedback control with inequality constraints, Int.Jnl.Contr., 37, 583-597.

[4] Clarke,D.W, Mohtadi,C, and Tuffs,P.S, 1987, Generalized predictive control, parts 1 and 2, Automatica, 23, 137-160.

[5] Cutler,C, 1983, Dynamic matrix control: an optimal multivariable control algorithm with constraints, Ph.D thesis, University of Houston.

[6] Cutler,C, and Hawkins,R, 1988, Application of a large predictive multivariable controller to a hydro-cracker second-stage reaction, Proc.ACC, 284-291.

[7] Gill, Murray and Wright, 1981, Practical Optimization, (New York: Academic Press).

[8] Golub,G.H, and Van Loan,C.F, 1983, Matrix Computations, (Baltimore,MD: Johns Hopkins UP).

[9] Griffiths,G.W, Dynamic simulation and control in the process industries, SAST Technical Bulletin No.5, SAST Ltd, UK, May 1992.

[10] Kohn,W, 1991, Declarative control architecture, Comm.ACM, 34,no.8, 64-79.

[11] Laban,M, On-line aircraft state and parameter estimation, TU Delft ???

[12] Landells,K, 1991, Application of MBPC to a distillation column at the Grangemouth oil refinery, Colloquium on Industrial Applications of MBPC, November, (London: IEE).

[13] Mayne,D.Q, and Michalska,H, 1990, Receding horizon control of nonlinear systems, IEEE Trans. Auto. Contr., 35, 814-824.

[14] Mohtadi,C, Shaw,S, and Fisher,D, 1992, Frequency response characteristics of MIMO GPC, Intl.Jnl.Control, 55, 887-900.

[15] Prett,D.M. and Garcia,C.E, 1988, Fundamental Process Control, (Boston: Butterworths).

[16] Richalet,J.A, et al, 1978, Model predictive heuristic control: applications to industrial processes, Automatica, 14, 413-???.

[17] Schinkel,W.M.M. and Sommelet,P, 1990, Model based predictive control in the petroleum industry, ??? Workshop on Computer Integrated Design of Controlled Industrial Systems, Paris, April 26-27.

[18] Yang,T.H, and Polak,E, 1991, Moving horizon control of nonlinear systems with input saturation, disturbances, and plant uncertainty, Memo. UCB/ERL M91/82, Electronics Research Lab., Univ. of California at Berkeley.

[19] Zafiriou,E, 1990, Robust model predictive control of processes with hard constraints, Computers Chem.Engng., 14, 359-371.