ECE 680 Fall 2011

Model-based Predictive Control (MPC)

by

Stanislaw H. Żak

1 Introduction

The model-based predictive control (MPC) methodology is also referred to as the moving horizon control or the receding horizon control. The idea behind this approach can be explained using an example of driving a car. The driver looks at the road ahead of him and taking into account the present state and the previous action predicts his action up to some distance ahead, which we refer to as the prediction horizon. Based on the prediction, the driver adjusts the driving direction. The MPC main idea is illustrated in Figure 1.

Figure 1: The driver predicts future travel direction based on the current state of the car and the current position of the steering wheel.

The MPC is constructed using control and optimization tools. The objective of this write-up is to introduce the reader to the linear MPC which refers to the family of MPC schemes in which linear models of the controlled objects are used in the control law synthesis.

\begin{matrix}c \\ ◯ \end{matrix}2009 by Stanislaw H. Żak

In the MPC approach, the current control action is computed on-line rather than using a pre-computed, off-line, control law.

A model predictive controller uses, at each sampling instant, the plant’s current input and output measurements, the plant’s current state, and the plant’s model to

• calculate, over a finite horizon, a future control sequence that optimizes a given performance index and satisfies constraints on the control action;

• use the first control in the sequence as the plant’s input.

The MPC strategy is illustrated in Figure 2, where N_{p} is the prediction horizon, u(t + k|t) is the predicted control action at t + k given u(t). Similarly, y(t + k|t) is the predicted output at t + k given y(t).

Figure 2: Controller action construction using model-based predictive control (MPC) approach.

2 Basic Structure of MPC

In Figure 3, we show a basic structure of an MPC-controlled plant, where we assume that the plant’s state is available to us.

Figure 3: State feedback model predictive controller.

3 From Continuous to Discrete Models

Our objective here is to present a method for constructing linear discrete-time models from given linear continuous-time models. The obtained discrete models will be used to perform computations to generate control commands.

We use a sample-and-hold device that transforms a continuous signal, f(t), into the staircase signal,

f(kh), kh ≤ t < (k + 1)h,

where h is the sampling period. The sample and zero-order hold (ZOH) operation is illustrated in Figure 4

Figure 4: Sample and zero-order hold (ZOH) element operating on a continuous function.

Suppose that we are given a continuous-time model,

ẋ(t) = Ax(t) + Bu(t), x_{0} = x(t_{0})
y(t) = Cx(t).

The solution to the state equation is

x(t) = e^{A(t−t_{0})}x(t_{0}) + ∫_{t_{0}}^{t} e^{A(t−τ)}Bu(τ)dτ.  (1)

We assume that the input to the system is generated by a sample-and-hold device and has the form,

u(t) = u(k), kh ≤ t < (k + 1)h.

Let t_{0} = kh and t = (k + 1)h and let us use shorthand notation,

x(kh) = x(k)

Then taking into account that u(k) is constant on the interval [kh, (k +1)h) we represent (1) as

x(k + 1) = e^{Ah}x(k) + ∫_{kh}^{(k+1)h} e^{A(kh+h−τ)}Bu(k)dτ
= e^{Ah}x(k) + ∫_{kh}^{(k+1)h} e^{A(kh+h−τ)}Bdτ u(k).  (2)

Consider now the second term on the right-hand side of the above equation. Let

η = kh + h − τ.

Then we can represent (2) as

x(k + 1) = e^{Ah}x(k) + ∫_{kh}^{(k+1)h} e^{A(kh+h−τ)}Bdτ u(k)
= e^{Ah}x(k) + ∫_{0}^{h} e^{Aη}Bdη u(k)
= Φx(k) + Γu(k),  (3)

where

Φ = e^{Ah} and Γ = ∫_{0}^{h} e^{Aη}Bdη.

The discrete output equation has the form

y(k) = Cx(k).

4 Simple Discrete-Time MPC

The development a discrete MPC controller in this section follows that of [9, Chapter 1]. We consider a discretized model of a dynamic system of the form,

x(k + 1) = Φx(k) + Γu(k)  (4)
y(k) = Cx(k),  (5)

where Φ ∈ ℝ^{n×n}, Γ ∈ ℝ^{n×m}, and C ∈ ℝ^{p×n}.

Applying the backward difference operator, Δx(k + 1) = x(k + 1) − x(k), to (4) gives

Δx(k + 1) = ΦΔx(k) + ΓΔu(k),  (6)

where Δu(k + 1) = u(k + 1) − u(k).

We now apply the backward difference operator to (5) to obtain

Δy(k + 1) = y(k + 1) − y(k)
= Cx(k + 1) − Cx(k)
= CΔx(k + 1).

Substituting into the above (6) yields

Δy(k + 1) = CΦΔx(k) + CΓΔu(k).

Hence,

y(k + 1) = y(k) + CΦΔx(k) + CΓΔu(k).  (7)

We combine (6) and (7) into one equation to obtain

[ \begin{matrix}Δx(k + 1) \\ y(k + 1) \end{matrix} ] = [ \begin{matrix}Φ & O \\ CΦ & I_{p} \end{matrix} ] [ \begin{matrix}Δx(k) \\ y(k) \end{matrix} ] + [ \begin{matrix}Γ \\ C & Γ \end{matrix} ] Δu(k).  (8)

We represent (5) as

y(k) = [ O I_{p} ] [ \begin{matrix}Δx(k) \\ y(k) \end{matrix} ] .  (9)

We now define the augmented state vector,

x_{a}(k) = [ \begin{matrix}Δx(k) \\ y(k) \end{matrix} ] .  (10)

Let

Φ_{a} = [ \begin{matrix}Φ & O \\ CΦ & I_{p} \end{matrix} ] , Γ_{a} = [ \begin{matrix}Γ \\ C & Γ \end{matrix} ] , and C_{a} = [ O I_{p} ] .  (11)

Using the above notation, we represent (8) and (9) in a compact format as

x_{a}(k + 1) = Φ_{a}x_{a}(k) + Γ_{a}Δu(k)  (12)
y(k) = C_{a}x_{a}(k),  (13)

where Φ_{a} ∈ ℝ^{(n+p)×(n+p)}, Γ_{a} ∈ ℝ^{(n+p)×m}, and C_{a} ∈ ℝ^{p×(n+p)}.

Suppose now that the state vector x_{a} at each sampling time, k, is available to us. Our control objective is to construct a control sequence,

Δu(k), Δu(k + 1), …, Δu(k + N_{p} − 1),  (14)

where N_{p} is the prediction horizon, such that a given cost function and constraints are satisfied. The above control sequence will result in a predicted sequence of the state vectors,

x_{a}(k + 1|k), x_{a}(k + 2|k), …, x_{a}(k + N_{p}|k),

which can then be used to compute predicted sequence of the plant’s outputs,

y(k + 1|k), y(k + 2|k), …, y(k + N_{p}|k).  (15)

Using the above information, we can compute the control sequence (14) and then apply u(k) to the plant to generate x(k + 1). We repeat the process again, using x(k + 1) as an initial condition to compute u(k + 1), until we reach the boundary of the control horizon, N_{c}.

We now present an approach to construct u(k) given x(k). Using the plant model parameters and the measurement of x_{a}(k) we evaluate the augmented states over the prediction horizon successively applying the recursion formula (12) to obtain,

x_{a}(k + 1|k) = Φ_{a}x_{a}(k) + Γ_{a}Δu(k)
x_{a}(k + 2|k) = Φ_{a}x_{a}(k + 1|k) + Γ_{a}Δu(k + 1)
= Φ_{a}^{2}x_{a}(k) + Φ_{a}Γ_{a}Δu(k) + Γ_{a}Δu(k + 1)
⋮
x_{a}(k + N_{p}|k) = Φ_{a}^{N_{p}}x_{a}(k) + Φ_{a}^{N_{p}−1}Γ_{a}Δu(k) + ⋯ + Γ_{a}Δu(k + N_{p} − 1)

We represent the above set of equations in the form,

[ x_{a}(k + 1|k) ; x_{a}(k + 2|k) ; ⋮ ; x_{a}(k + N_{p}|k) ] = [ Φ_{a} ; Φ_{a}^{2} ; ⋮ ; Φ_{a}^{N_{p}} ] x_{a}(k) + [ \begin{matrix}Γ & a \\ Φ_{a}Γ_{a} & Γ_{a} \\ ⋮ & & ... \\ Φ_{a}^{N_{p}−1}Γ_{a} & & ⋯ & Γ_{a} \end{matrix} ] [ Δu(k) ; Δu(k + 1) ; ⋮ ; Δu(k + N_{p} − 1) ] .

(16)

We wish to design a controller that would force the plant output, y, to track a given reference signal, r. Using (13) and (16), we compute the sequence of predicted outputs (15),

[ y(k + 1|k) ; y(k + 2|k) ; ⋮ ; y(k + N_{p}|k) ] = [ C_{a}x_{a}(k + 1|k) ; C_{a}x_{a}(k + 2|k) ; ⋮ ; C_{a}x_{a}(k + N_{p}|k) ]

= [ C_{a}Φ_{a} ; C_{a}Φ_{a}^{2} ; ⋮ ; C_{a}Φ_{a}^{N_{p}} ] x_{a}(k)  (17)

+ [ \begin{matrix}C & a & Γ & a \\ C_{a}Φ_{a}Γ_{a} & C_{a}Γ_{a} \\ ⋮ & & ... \\ C_{a}Φ_{a}^{N_{p}−1}Γ_{a} & & ⋯ & C_{a}Γ_{a} \end{matrix} ] [ Δu(k) ; Δu(k + 1) ; ⋮ ; Δu(k + N_{p} − 1) ] .  (18)

We write the above compactly as

Y = W x_{a}(k) + ZΔU,  (19)

where

Y = [ y(k + 1|k) ; y(k + 2|k) ; ⋮ ; y(k + N_{p}|k) ] , ΔU = [ Δu(k) ; Δu(k + 1) ; ⋮ ; Δu(k + N_{p} − 1) ] ,

and

W = [ C_{a}Φ_{a} ; C_{a}Φ_{a}^{2} ; ⋮ ; C_{a}Φ_{a}^{N_{p}} ] , and Z = [ \begin{matrix}C & a & Γ & a \\ C_{a}Φ_{a}Γ_{a} & C_{a}Γ_{a} \\ ⋮ & & ... \\ C_{a}Φ_{a}^{N_{p}−1}Γ_{a} & & ⋯ & C_{a}Γ_{a} \end{matrix} ] .

Suppose now that we wish to construct a control sequence, Δu(k), …, Δu(k + N_{p} − 1), that would minimize the cost function

J(ΔU) = \frac{1}{2} (r_{p} − Y )^{⊤} Q (r_{p} − Y ) + \frac{1}{2}ΔU^{⊤}RΔU,  (20)

where Q = Q^{⊤} > 0 and R = R^{⊤} > 0 are real symmetric positive semi-definite and positivedefinite weight matrices, respectively. The multiplying scalar, 1/2, is just to make subsequent manipulations cleaner. Finally, the vector r_{p} consists of the values of the command signal at sampling times, k, k + 1, …, k + N_{p} − 1. The selection of the weight matrices, Q and R reflects our control objective to keep the tracking error ∥r_{p} − Y ∥ “small” using the control actions that are “not too large.”

We first apply the first-order necessary condition (FONC) test to J(ΔU),

\frac{∂J}{∂ΔU} = 0^{⊤}.

Then, we solve the above equation for ΔU = ΔU^{∗}, where

\frac{∂J}{∂ΔU} = − (r_{p} − W x_{a} − ZΔU)^{⊤} QZ + ΔU^{⊤}R
= 0^{⊤}.

Performing simple manipulations yields

−r_{p}^{⊤}QZ + x_{a}^{⊤}W^{⊤}QZ + ΔU^{⊤}Z^{⊤}QZ + ΔU^{⊤}R = 0^{⊤}.

Applying the transposition operation to both sides of the above equation and rearranging terms, we obtain

(R + Z^{⊤}QZ) ΔU = Z^{⊤}Q (r_{p} − W x_{a}) .

Note that the matrix (R + Z^{⊤}QZ) is invertible, and in fact, positive definite because R = R^{⊤} > 0 and Z^{⊤}QZ is also symmetric and at least positive semi-definite. Hence, ΔU that satisfies the FONC is

ΔU^{∗} = (R + Z^{⊤}QZ)^{−1} Z^{⊤}Q (r_{p} − W x_{a}) .  (21)

Now, applying the second derivative test to J(ΔU), which we refer to as the second-order sufficiency condition (SONC), we obtain

\frac{∂^{2}J}{∂ΔU^{2}} = R + Z^{⊤}QZ
> 0,

Figure 5: Discrete-time MPC for linear time-invariant systems.

which implies that ΔU^{∗} is a strict minimizer of J.

Using (21), we compute Δu(k),

Δu(k) = \overbrace{[ I_{m} O ⋯ O ]}^{N_{p} block matrices} (R + Z^{⊤}QZ)^{−1} Z^{⊤}Q (r_{p} − W x_{a})
= Krr_{p} − KxΔx(k) − Kyy(k),  (22)

where

Kr = [ I_{m} O ⋯ O ] (R + Z^{⊤}QZ)^{−1} Z^{⊤}Q, Kx = KrW [ \begin{matrix}I & n \\ O \end{matrix} ] ,

and

Ky = KrW [ \begin{matrix}O \\ I & p \end{matrix} ] .

An implementation of the above controller, using a discrete-time integrator, is shown in Figure 5.

5 MPC With Constraints

An attractive feature of the model-based predictive control approach is that a control engineer can incorporate different types of constraints on the control action. We consider three types of such constraints.

5.1 Constraints on the Rate of Change of the Control Action

Hard constraints on the rate of change of the control signal can be expressed as

Δu_{i}^{min} ≤ Δu_{i}(k) ≤ Δu_{i}^{max}, i = 1, 2, …, m  (23)

Let

Δu^{min} = [ Δu_{1}^{min} ⋯ Δu_{m}^{min} ]^{⊤} and Δu^{max} = [ Δu_{1}^{max} ⋯ Δu_{m}^{max} ]^{⊤} .

Then, we can express (23) as

Δu^{min} ≤ Δu(k) ≤ Δu^{max}.  (24)

The above can be equivalently represented as

[ \begin{matrix}−I_{m} \\ I & m \end{matrix} ] Δu(k) ≤ [ \begin{matrix}−Δu^{min} \\ Δu^{max} \end{matrix} ] .  (25)

Using the above approach we can represent constraints on the rate of change of the control action over the whole prediction horizon, N_{p}, in terms of ΔU(k), by augmenting the above inequality to incorporate constraints for the remaining sampling times. That is, if the above constraints are imposed on the rate of change of the control action for all sampling times within the prediction horizon, then this can be expressed in terms of ΔU(k) as

[ \begin{matrix}−I_{m} & O & ⋯ & O & O \\ I_{m} & O & ⋯ & O & O \\ O & −I_{m} & ⋯ & O & O \\ O & I_{m} & ⋯ & O & O \\ ⋮ & & & & ⋮ \\ O & O & ⋯ & −I_{m} & O \\ O & O & ⋯ & I_{m} & O \\ O & O & ⋯ & O & −I_{m} \\ O & O & ⋯ & O & I_{m} \end{matrix} ] [ Δu(k) ; Δu(k + 1) ; ⋮ ; Δu(k + N_{p} − 2) ; Δu(k + N_{p} − 1) ] ≤ [ −Δu^{min} Δu^{max} −Δu^{min} Δu^{max} ⋮ −Δu^{min} Δu^{max} −Δu^{min} Δu^{max} ] .  (26)

On the other hand, if the constraints on the rate of change of the control are imposed only on the first component of ΔU(k), then we express this in terms of ΔU(k) as

[ \begin{matrix}−I_{m} & O & ⋯ & O & O \\ I_{m} & O & ⋯ & O & O \end{matrix} ] ΔU(k) ≤ [ \begin{matrix}−Δu^{min} \\ Δu^{max} \end{matrix} ] .

5.2 Constraints on the Control Action Magnitude

Hard constraints on the control action magnitude at the time sampling k have the form

u_{i}^{min} ≤ u_{i}(k) ≤ u_{i}^{max}, i = 1, 2 …, m  (27)

We now express constraints on the control action magnitude over the whole prediction horizon in term of ΔU(k). To accomplish our goal, we first note that

u(k) = u(k − 1) + Δu(k)
= u(k − 1) + [ I_{m} O ⋯ O ] ΔU(k),

where, recall that ΔU(k) = [ Δu(k) Δu(k + 1) ⋯ Δu(k + N_{p} − 1) ]^{⊤}.

Similarly,

u(k + 1) = u(k) + Δu(k + 1)
= u(k − 1) + [ I_{m} O ⋯ O ] ΔU(k) + Δu(k + 1)
= u(k − 1) + [ I_{m} I_{m} ⋯ O ] ΔU(k).

Continuing in this manner, we obtain

[ u(k) ; u(k + 1) ; ⋮ ; u(k + N_{p} − 1) ] = [ I_{m} ; I_{m} ; ⋮ ; I_{m} ] u(k − 1) + [ \begin{matrix}I_{m} & O & ⋯ & O \\ I_{m} & I_{m} & ⋯ & O \\ ⋮ & & ... & ⋮ \\ I_{m} & I_{m} & ⋯ & I_{m} \end{matrix} ] [ Δu(k) ; Δu(k + 1) ; ⋮ ; Δu(k + N_{p} − 1) ] . (28)

Let

E = [ I_{m} ; I_{m} ; ⋮ ; I_{m} ] and H = [ \begin{matrix}I_{m} & O & ⋯ & O \\ I_{m} & I_{m} & ⋯ & O \\ ⋮ & & ... & ⋮ \\ I_{m} & I_{m} & ⋯ & I_{m} \end{matrix} ] .

Then, we can represent (28) as

U(k) = Eu(k − 1) + HΔU(k).  (29)

Suppose now that we are faced with constructing a control action subject to the following constraints,

U^{min} ≤ U(k) ≤ U^{max}.

The above constraints can be equivalently represented as

[ \begin{matrix}−U(k) \\ U(k) \end{matrix} ] ≤ [ \begin{matrix}−U^{min} \\ U & max \end{matrix} ] ,

Taking into account (29), we write the above as

[ \begin{matrix}− (Eu(k − 1) + HΔU(k)) \\ Eu(k − 1) + HΔU(k) \end{matrix} ] ≤ [ \begin{matrix}−U^{min} \\ U & max \end{matrix} ] .  (30)

The above, in turn, can be represented as

[ \begin{matrix}−H \\ H \end{matrix} ] ΔU(k) ≤ [ \begin{matrix}−U^{min} + Eu(k − 1) \\ U^{max} − Eu(k − 1) \end{matrix} ] .  (31)

In a special case, when a control designer elects to impose constraints only on the first component of ΔU(k), that is, on Δu(k) only, then this scenario is expressed in terms of ΔU(k) as

[ \begin{matrix}−I_{m} & O & ⋯ & O \\ I_{m} & O & ⋯ & O \end{matrix} ] ΔU(k) ≤ [ \begin{matrix}−U^{min} + Eu(k − 1) \\ U^{max} − Eu(k − 1) \end{matrix} ] .  (32)

5.3 Constraints on the Plant Output

Recall from (19) that the predicted plant output is, Y (k) = W x_{a}(k) + ZΔU(k). Suppose now that the following constraints are imposed on the predicted plant’s output,

Y^{min} ≤ Y (k) ≤ Y^{max}.

We represent the above as

[ \begin{matrix}−Y (k) \\ Y (k) \end{matrix} ] ≤ [ \begin{matrix}−Y^{min} \\ Y & max \end{matrix} ] .

Taking into the account the expression for Y given by (19), we represent the above as

[ \begin{matrix}−Z \\ Z \end{matrix} ] ΔU(k) ≤ [ \begin{matrix}−Y^{min} + W x_{a}(k) \\ Y^{max} − W x_{a}(k) \end{matrix} ] .  (33)

The above discussion clearly demonstrates that we need an effective method of minimizing a function of many variables, J(ΔU), subject to inequality constraints such as given by (25), (30), and (33), or their combination. In the following, we present a powerful method for solving such optimization problems. This is the method that we shall use to implement our MPCs.

6 An Optimizer for Solving Constrained Optimization Problems

If follows from the discussion in the previous section that at each sampling time the MPC calls for a solution to a constrained optimization problem of the form,

minimize J(ΔU)
subject to g(ΔU) ≤ 0,

where g(ΔU) ≤ 0 contains inequality constrains given by (25), (30), and (33), or their combination. In this section, we present an iterative methods for solving the above optimization problems. To proceed, we first present the descent gradient method, followed by the Newton’s method, for solving unconstrained optimization problems of the form,

minimize J(ΔU)

6.1 Gradient Descent Method

For the sake of simplicity, we denote the argument of a function of many variables as x, where x ∈ ℝ^{N}, and the function that we will be minimizing will be denoted as f, where f : ℝ^{N} → ℝ.

The method of the gradient descent is based on the following property of the gradient of a differentiable function, f, on ℝ^{N}:

Theorem 1 At a given point x^{(0)}, the vector

v = −∇f (x^{(0)})

points in the direction of most rapid decrease of f and the rate of increase of f at x^{(0)} in the direction v is − ‖∇f (x^{(0)})‖, equivalently, the rate of decrease of f at x^{(0)} in the direction v is ‖∇f (x^{(0)})‖.

Thus, if we wish to minimize a differentiable function, f, then moving in the direction of the negative gradient is a good direction. The gradient descent algorithm rests on the above observation and has the form,

x^{(k+1)} = x^{(k)} − α∇f (x^{(k)}) ,

where α > 0 is a step size.

6.2 Second-Order Necessary Conditions for a Minimum

Suppose we are given a function f of one variable x. Recall a well-known theorem referred to as the second-order Taylor’s formula or the extended law of mean [8, p. 1]:

Theorem 2 Suppose that f(x), f′(x), f′′(x) exist on the closed interval [a, b] = {x ∈ ℝ : a ≤ x ≤ b}. If x^{∗}, x are any two different points of [a, b], then there exists a point z strictly between x^{∗} and x such that

f(x) = f (x^{∗}) + f′ (x^{∗}) (x − x^{∗}) + \frac{f′′(z)}{2} (x − x^{∗})^{2} .

Using the above formula we observe that if

• f′ (x^{∗}) = 0, and

• f′′ (x^{∗}) > 0,

then

f(x) = f (x^{∗}) + a positive number

for all x “close” to x^{∗}. Indeed, if f′′(x) is continuous at x^{∗} and f′′ (x^{∗}) > 0, then f′′ (x) > 0 for all x in some neighborhood of x^{∗}. Therefore,

f(x) > f (x^{∗}) for all x close to x^{∗},

which means that x^{∗} is a strict local minimizer of f.

We now extend the above result for functions of many variables. To proceed, we need the second-order Taylor’s formula for such functions that can be found in [8, p. 11].

Theorem 3 Suppose that x^{∗}, x are points in ℝ^{N} and that f is a function of N variables with continuous first and second partial derivatives on some open set containing the line segment

[x^{∗}, x] = {w ∈ ℝ^{N} : w = x^{∗} + t (x − x^{∗}) ; 0 ≤ t ≤ 1}

joining x^{∗} and x. Then there exists a z ∈ [x^{∗}, x] such that

f(x) = f (x^{∗}) + ∇f (x^{∗})^{⊤} (x − x^{∗}) + \frac{1}{2} (x − x^{∗})^{⊤} F (z) (x − x^{∗}) ,  (34)

where F (·) is the Hessian of f, that is, the second derivative of f.

Hence if

• x^{∗} is a critical point, that is, ∇f (x^{∗}) = 0, and

• F (x^{∗}) > 0,

then using (34), we conclude that

f(x) = f (x^{∗}) + 0 + a positive number

for all x in a neighborhood of x^{∗}. Therefore for all x ≠ x^{∗} in some neighborhood of x^{∗}, we have

f(x) > f (x^{∗}) ,

which implies that x^{∗} is a strict local minimizer.

6.3 Newton’s Method

The idea behind Newton’s method for function minimization is to minimize the quadratic approximation rather than the function itself as illustrated in Figure 6. Newton’s method

Figure 6: Newton’s method minimizes the quadratic approximation of the objective function that utilizes first and second derivatives of the optimized function.

seeks a critical point, x^{∗}, of a given function. If at this critical point we have F (x^{∗}) > 0, then x^{∗} is a strict local minimizer of f.

We can obtain a quadratic approximation q of f at x^{∗} from the second-order Taylor series expansion of f about x^{∗},

q(x) = f (x^{∗}) + ∇f (x^{∗})^{⊤} (x − x^{∗}) + \frac{1}{2} (x − x^{∗})^{⊤} F (x^{∗}) (x − x^{∗}) .

Note that

q (x^{∗}) = f (x^{∗}) , ∇q (x^{∗}) = ∇f (x^{∗}) ,

as well as their Hessians, that is, their second derivatives evaluated at x^{∗} are equal. A critical point of q can be obtained by solving the algebraic equation,

∇q(x) = 0,

that is, by solving the equation

∇q(x) = ∇f (x^{∗}) + F (x^{∗}) (x − x^{∗}) = 0.

Suppose now that we have a quadratic approximation of f at a point x^{(k)}, that is,

q(x) = f (x^{(k)}) + ∇f (x^{(k)})^{⊤} (x − x^{(k)}) + \frac{1}{2} (x − x^{(k)})^{⊤} F (x^{(k)}) (x − x^{(k)}) .

We assume that det F (x^{(k)}) ≠ 0. Denoting the solution to the above equation as x^{(k+1)}, we obtain

x^{(k+1)} = x^{(k)} − F (x^{(k)})^{−1} ∇f (x^{(k)})

The above is known as the Newton’s method for minimizing a function of many variables f. Note that x^{(k+1)} is a critical point of the quadratic function q that approximates f at x^{(k)}.

A computationally efficient representation of Newton’s algorithm has the form,

x^{(k+1)} = x^{(k)} − Δx^{(k)},

where Δx^{(k)} is obtained by solving the equation,

F (x^{(k)}) Δx^{(k)} = ∇f (x^{(k)}) .

7 Minimization Subject to Equality Constraints

We now discuss the problem of finding a point x ∈ ℝ^{N} that minimizes f(x) subject to equality constraints,

h_{1}(x) = 0 ; h_{2}(x) = 0 ; ⋮ ; h_{M}(x) = 0 }

where M ≤ N. We write the above equality constraints in a compact form as

h(x) = 0,

where h : ℝ^{N} → ℝ^{M}. We refer to the set of points satisfying the above constraints as the surface.

We now introduce the notion of the tangent plane to the surface S = {x : h(x) = 0} at a point x^{∗} ∈ S . First, we define a curve on the surface S as a family of points x(t) ∈ S continuously parameterized by t for t ∈ [a, b]. The curve is differentiable if ẋ(t) = dx(t)/dt exists. A curve is said to pass through the point x^{∗} ∈ S if x^{∗} = x (t^{∗}) for some t^{∗} ∈ [a, b]. The tangent plane to the surface S = {x : h(x) = 0} at x^{∗} is the collection of the derivatives at x^{∗} of all differentiable curves on S that pass through x^{∗}.

To proceed, we need the notion of a regular point of the constraints. We say that x^{∗} satisfying the constraints, that is, h (x^{∗}) = 0, is a regular point of the constraints if the gradient vectors,

∇h_{1} (x^{∗}) , …, ∇h_{M} (x^{∗})

are linearly independent.

One can show, see, for example [4, p. 298], that the tangent space at a regular point x^{∗}, denoted T (x^{∗}), to the surface {x : h (x) = 0} at the regular point x^{∗} is

T (x^{∗}) = {y : [ ∇h_{1} (x^{∗})^{⊤} ; ⋮ ; ∇h_{M} (x^{∗})^{⊤} ] y = 0} .  (35)

With the above notions in place, we are ready to state and prove the following result which is known as the first-order necessary condition (FONC) for function minimization subject to equality constraints.

Theorem 4 Let x^{∗} be a local minimizer (or maximizer) of f subject to the constraints h(x) = 0 and let x^{∗} be a regular point of the constraints. Then there exists a vector λ^{∗} such that

∇f (x^{∗}) + [ ∇h_{1} (x^{∗}) ⋯ ∇h_{M} (x^{∗}) ] λ^{∗} = 0.

Proof Let x(t) be a differentiable curve passing through x^{∗} on the surface S = {x : h(x) = 0} such that ẋ (t^{∗}) = y where t^{∗} ∈ [a, b]. Note that the tangent plane is given by (35). Because x^{∗} is a local minimizer of f on S, we have

\frac{d}{dt}f(x(t))|_{t=t^{∗}} = 0.

Applying the chain rule to the above gives

∇f (x^{∗})^{⊤} y = 0.

Thus ∇f (x^{∗}) is orthogonal to the tangent space T (x^{∗}). That is, ∇f (x^{∗}) is a linear combination of the gradients ∇h_{1} (x^{∗}) , …, ∇h_{M} (x^{∗}). This fact can be expressed as

∇f (x^{∗}) + [ ∇h_{1} (x^{∗}) ⋯ ∇h_{M} (x^{∗}) ] λ^{∗} = 0

for some constant vector λ^{∗} ∈ ℝ^{M}.

□

The vector λ^{∗} is called the vector of Lagrange multipliers.

We now introduce the Lagrangian associated with the constrained optimization problem,

l(x, λ) = f(x) + λ^{⊤}h(x).

Then the FONC can be expressed as

∇xl(x, λ) = 0
∇λl(x, λ) = 0.

Note that the second of the above condition is equivalent to h(x) = 0, that is,

∇λl(x, λ) = h(x) = 0.

Equivalently the FONC can be written as

∇xl(x, λ) = 0
h(x) = 0.

We now apply Newton’s method to solve the above system of equations iteratively,

[ \begin{matrix}x & (k+1) \\ λ & (k+1) \end{matrix} ] = [ \begin{matrix}x & (k) \\ λ & (k) \end{matrix} ] + [ \begin{matrix}d & (k) \\ y & (k) \end{matrix} ] ,

where d^{(k)} and y^{(k)} are obtained by solving the matrix equation,

[ \begin{matrix}L (x^{(k)}, λ^{(k)}) & Dh (x^{(k)})^{⊤} \\ Dh (x^{(k)}) & O \end{matrix} ] [ \begin{matrix}d & (k) \\ y & (k) \end{matrix} ] = [ \begin{matrix}−∇xl(x, λ) \\ −h (x^{(k)}) \end{matrix} ] ,

where L (x, λ) is the Hessian of l(x, λ) with respect to x, and Dh(x) is the Jacobian matrix of h(x), that is,

Dh(x) = [ ∇h_{1} (x)^{⊤} ; ⋮ ; ∇h_{M} (x)^{⊤} ] .

The above algorithm is also referred to as in the literature as sequential quadratic programming (SQP); see, for example [7, Section 15.5].

8 A Lagrangian Algorithm for Equality Constraints

The first-order Lagrangian algorithm for the optimization problem involving minimizing f subject to the equality constraints, h(x) = 0, has the form,

x^{(k+1)} = x^{(k)} − α_{k} (∇f (x^{(k)}) + Dh (x^{(k)})^{⊤} λ^{(k)})
λ^{(k+1)} = λ^{(k)} + β_{k}h (x^{(k)}) ,

where α_{k} and β_{k} are positive constants. Note that the update for x^{(k)} is a descent gradient for minimizing the Lagrangian with respect to x, while the update for λ^{(k)} is a gradient ascent for maximizing the Lagrangian with respect to λ. For an analysis of the above algorithm, we recommend [2, pp. 522–524]

9 Minimization Subject to Inequality Constraints

We now discuss the problem of finding a point x ∈ ℝ^{N} that minimizes f(x) subject to inequality constraints,

g_{1}(x) ≤ 0 ; g_{2}(x) ≤ 0 ; ⋮ ; g_{P} (x) ≤ 0 }

We write the above inequality constraints in a compact form as

g(x) ≤ 0,

where h : ℝ^{N} → ℝ^{P}.

Let x^{∗} be a point satisfying the constraints, that is, g (x^{∗}) ≤ 0, and let J be the set of indices j for which g_{j} (x^{∗}) = 0. Then x^{∗} is said to be a regular point of the constraints g(x) ≤ 0 if the gradient vectors,

∇g_{j} (x^{∗}) , j ∈ J,

are linearly independent.

We say that a constraint g_{j}(x) ≤ 0 is active at x^{∗} if g_{j} (x^{∗}) = 0. Thus the index set J, defined above, contains indices of active constraints.

Theorem 5 Let x^{∗} be a regular point and a local minimizer of f subject to g(x) ≤ 0. Then there exists a vector μ^{∗} ∈ ℝ^{P} such that

1. μ^{∗} ≥ 0,

2. ∇f (x^{∗}) + [ ∇g_{1} (x^{∗}) ⋯ ∇g_{P} (x^{∗}) ] μ^{∗} = 0,

3. μ^{∗⊤}g (x^{∗}) = 0.

Proof First note that because x^{∗} is a relative minimizer over the constraint set {x : g(x) ≤ 0}, it is also a minimizer over a subset of the constraint set obtained by setting the active constraints to zero. Therefore, for the resulting equality constrained problem, by Theorem 4, we have

∇f (x^{∗}) + [ ∇g_{1} (x^{∗}) ⋯ ∇g_{P} (x^{∗}) ] μ^{∗} = 0,  (36)

where μ_{j}^{∗} = 0 if g_{j} (x^{∗}) < 0. This means that μ^{∗⊤}g (x^{∗}) = 0. Thus μ_{i} may be non-zero only if the corresponding constraint is active, that is, g_{i} (x^{∗}) = 0.

It thus remains to show that μ^{∗} ≥ 0. We prove this by contraposition. We suppose that for some k ∈ J, we have μ_{k}^{∗} < 0. Consider next a surface formed by all other active constraints, that is, the surface

{x : g_{j}(x) = 0, j ∈ J, j ≠ k}.  (37)

Also consider a tangent space to the above surface at x^{∗}. By assumption x^{∗} is a regular point of the active constraints. Therefore, there exists a vector y such that

∇g_{k} (y)^{⊤} y < 0.

Let now x(t), t ∈ [−a, a], a > 0, be a curve on the surface (37) such that

ẋ(0) = y.

Note that for small t, the curve x(t) is feasible. We apply the transposition operator to the both sides of (36) to obtain

∇f (x^{∗})^{⊤} + ∑_{i=1}^{P} μ_{i}^{∗}∇g_{i} (x^{∗})^{⊤} = 0^{⊤}.

Post-multiplying the above by y an taking into account the fact that y belongs to the tangent space to the surface (37) gives

∇f (x^{∗})^{⊤} y = −μ_{k}^{∗}∇g_{k} (x^{∗})^{⊤} y
< 0.

Suppose, without loss of generality, that ∥y∥_{2} = 1. Then, we have that

\frac{df(x(t))}{dt} = ∇f (x^{∗})^{⊤} y
< 0,

that is, the rate of increase of f at x^{∗} in the direction y is negative. This would mean that we could decrease the value of f moving just slightly away from x^{∗} along y while, at the same time, preserving feasibility. But this contradicts the minimality of f (x^{∗}). In summary, if x^{∗} is a relative minimizer then we also have the components of μ^{∗} all non-negative.

□

The vector μ^{∗} is called the vector of the Karush-Kuhn-Tucker (KKT) multipliers.

10 A Lagrangian Algorithm for Inequality Constraints

We now present a first-order Lagrangian algorithm for the optimization problem involving inequality constraints,

minimize f(x)
subject to g(x) ≤ 0,

where g : ℝ^{N} → g : ℝ^{P}. The Lagrangian function is

l(x, μ) = f(x) + μ^{⊤}g(x).

The first-order Lagrangian algorithm for the above optimization problem involving minimizing f subject to the inequality constraints, g(x) ≤ 0, has the form,

x^{(k+1)} = x^{(k)} − α_{k} (∇f (x^{(k)}) + Dg (x^{(k)})^{⊤} μ^{(k)})
μ^{(k+1)} = [μ^{(k)} + β_{k}g (x^{(k)})]_{+} ,

where the operation [·]_{+} = max(·, 0) is applied component-wise.

For an analysis of the above algorithm, we recommend [2, pp. 524–528]

Example 1 We illustrate the design of the discrete MPC subject to the constraints on the control input on a plant modeled by the following continuous state-space model,

\begin{matrix}ẋ & = & [ \begin{matrix}−0.1 & −3.0 \\ 1 & 0 \end{matrix} ] x + [ \begin{matrix}1 \\ 0 \end{matrix} ] u \\ y & = & [ 0 10 ] x. \end{matrix} \}  (38)

The system with the transfer function as the above system was used by Wang [9, p. 70] to test her discrete MPC designs. We use different optimizer than Wang. The reader can compare the results of two different MPC synthesis methods. The MPC controller’s objective is to force the plant’s output to track the unit step.

We proceed to generate the discrete state-space model employing the MATLAB’s command c2dm. We use the sampling time h = 0.1. The command input is the unit step. The weight matrices in the cost function (20) are R = 0.01I_{3} and Q = I_{3}. The prediction horizon is N_{p} = 3. Our objective is to design a discrete-MPC subject to the constraints on the control of the form,

−0.3 ≤ u(k) ≤ 0.5

First, we consider the case when the constraints are imposed only on the first component of U(k). The constraints take the form of (32). Thus, in our example, the constraints can be written as,

[ \begin{matrix}−1 & 0 & 0 \\ 1 & 0 & 0 \end{matrix} ] [ \begin{matrix}Δu(k) \\ Δu(k + 1) \\ Δu(k + 2) \end{matrix} ] ≤ [ \begin{matrix}0.3 + u(k − 1) \\ 0.5 − u(k − 1) \end{matrix} ] .

We employed the first-order Lagrangian algorithm for inequality constraints presented in Section 10, where α_{k} = β_{k} = 0.005. We set zero initial conditions. Plots of the control effort and the plant’s output are shown in Figure 7. We wish to emphasize that the simulations were performed for the discretized plant. Thus, the control and output were evaluated at the sampling instances and marked on the plots using circles. For the purpose of visualization the control action was presented in the form of a staircase function while the plant output values were connected by straight lines.

Next we tested the discrete MPC where the constraints, −0.3 ≤ u(k) ≤ 0.5, were imposed Figure 7: Plots of control effort and the plant’s output versus time of the closed-loop system of Example 1 for the case when the constraint was imposed only on the first component of U(k).

on all the elements of U(k). We use (31) to express the constraints, where

[ \begin{matrix}−1 & 0 & 0 \\ −1 & −1 & 0 \\ −1 & −1 & −1 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{matrix} ] [ \begin{matrix}Δu(k) \\ Δu(k + 1) \\ Δu(k + 2) \end{matrix} ] ≤ [ \begin{matrix}0.3 + u(k − 1) \\ 0.3 + u(k − 1) \\ 0.3 + u(k − 1) \\ 0.5 − u(k − 1) \\ 0.5 − u(k − 1) \\ 0.5 − u(k − 1) \end{matrix} ] .

Example 2 In this example, we apply the discrete MPC controller from Example 1 to the continuous plant modeled by (38). The MPC controller’s objective is to force the plant’s output to track the unit step. The zero-order hold is applied to the discrete MPC controller’s output sequence resulting in a piece-wise constant input to the continuous plant model. In Figure 9, we show plots of the plant output as well as the control effort versus time. As in the previous example, the prediction horizon is N_{p} = 3 and the constraints on the control are

−0.3 ≤ u(k) ≤ 0.5

Figure 8: Plots of control effort and the plant’s output versus time of the closed-loop system of Example 1 for the case when the constraints were imposed on all the components of U(k).

11 Observer-Based MPC

If the plant state is not accessible, we need to use an observer in the controller implementation as shown in Figure 10.

12 Notes

The key reason for huge popularity of the MPC approach is its ability to systematically take into account constraints thus allowing processes to operate at the limits of achievable performance [3]. For some impressive industrial applications of MPCs see [3, 1]. Nonlinear model predictive controllers are presented in [5, 6].

Figure 9: Plots of control effort and the plant’s output versus time of the closed-loop system of Example 2 for the case when the constraints were imposed only on the first component of U(k).

References

[1] E. F. Camacho and C. Bordons. Model Predictive Control. Springer, London, second edition, 2004.

[2] E. K. P. Chong and S. H. Żak. An Introduction to Optimization. Wiley-Interscience, John Wiley & Sons, Inc., New York, third edition, 2008.

[3] B. Kouvaritakis and M. Cannon, editors. Nonlinear Predictive Control: Theory and Practice. The Institution of Electrical Engineers, London, 2001.

[4] D. G. Luenberger. Linear and Nonlinear Programming. Addison-Wesley Publishing, Reading, Massachusetts, second edition, 1984.

[5] D. Q. Mayne and H. Michalska. Receding horizon control of nonlinear systems. IEEE Transactions on Automatic Control, 35(7):814–824, July 1990.

[6] H. Michalska and D. Q. Mayne. Moving horizon observers and observer-based control. IEEE Transactions on Automatic Control, 40(6):995–1006, June 1995.

Figure 10: MPC implementation when the plant’s state is not available.

[7] S. G. Nash and A. Sofer. Linear and Nonlinear Programming. McGraw-Hill, New York, 1996.

[8] A. L. Peressini, F. E. Sullivan, and J. J. Uhl, Jr. The Mathematics of Nonlinear Programming. Springer-Verlag, New York, 1988.

[9] L. Wang. Model Predictive Control System Design and Implementation Using MATLAB. Springer, London, 2009.