Optimization Problem

Overview

Model predictive control solves an optimization problem – specifically, a quadratic program (QP) – at each control interval. The solution determines the manipulated variables (MVs) to be used in the plant until the next control interval.

This QP problem includes the following features:

The objective, or "cost", function — A scalar, nonnegative measure of controller performance to be minimized.
Constraints — Conditions the solution must satisfy, such as physical bounds on MVs and plant output variables.
Decision — The MV adjustments that minimize the cost function while satisfying the constraints.

The following sections describe these features in more detail.

Standard Cost Function

The standard cost function is the sum of four terms, each focusing on a particular aspect of controller performance, as follows:

$J (z_{k}) = J_{y} (z_{k}) + J_{u} (z_{k}) + J_{Δ u} (z_{k}) + J_{ε} (z_{k}) .$

Here, z_k is the QP decision. As described below, each term includes weights that help you balance competing objectives. While the MPC controller provides default weights, you will usually need to adjust them to tune the controller for your application.

Output Reference Tracking

In most applications, the controller must keep selected plant outputs at or near specified reference values. An MPC controller uses the following scalar performance measure for output reference tracking:

$J_{y} (z_{k}) = \sum_{j = 1}^{n_{y}} \sum_{i = 1}^{p} {\frac{w_{i, j}^{y}}{s_{j}^{y}} [r_{j} (k + i | k) - y_{j} (k + i | k)]}^{2} .$

Here,

k — Current control interval.
p — Prediction horizon (number of intervals).
n_y — Number of plant output variables.
z_k — QP decision, given by:
$z_{k}^{T} = [\begin{matrix} u {(k | k)}^{T} & u {(k + 1 | k)}^{T} & \dots & u {(k + p - 1 | k)}^{T} & ε_{k} \end{matrix}] .$
y_j(k+i|k) — Predicted value of jth plant output at ith prediction horizon step, in engineering units.
r_j(k+i|k) — Reference value for jth plant output at ith prediction horizon step, in engineering units.
$s_{j}^{y}$ — Scale factor for jth plant output, in engineering units.
$w_{i, j}^{y}$ — Tuning weight for jth plant output at ith prediction horizon step (dimensionless).

The values n_y, p, $s_{j}^{y}$ , and $w_{i, j}^{y}$ are constant controller specifications. The controller receives reference values, r_j(k+i|k), for the entire prediction horizon. The controller uses the state observer to predict the plant outputs, y_j(k+i|k), which depend on manipulated variable adjustments (z_k), measured disturbances (MD), and state estimates. At interval k, the controller state estimates and MD values are available. Therefore, J_y is a function of z_k only.

Manipulated Variable Tracking

In some applications, such as when there are more manipulated variables than plant outputs, the controller must keep selected manipulated variables (MVs) at or near specified target values. An MPC controller uses the following scalar performance measure for manipulated variable tracking:

$J_{u} (z_{k}) = \sum_{j = 1}^{n_{u}} \sum_{i = 0}^{p - 1} {\frac{w_{i, j}^{u}}{s_{j}^{u}} [u_{j} (k + i | k) - u_{j, t a r g e t} (k + i | k)]}^{2} .$

Here,

k — Current control interval.
p — Prediction horizon (number of intervals).
n_u — Number of manipulated variables.
z_k — QP decision, given by:
$z_{k}^{T} = [\begin{matrix} u {(k | k)}^{T} & u {(k + 1 | k)}^{T} & \dots & u {(k + p - 1 | k)}^{T} & ε_{k} \end{matrix}] .$
u_j,target(k+i|k) — Target value for jth MV at ith prediction horizon step, in engineering units.
$s_{j}^{u}$ — Scale factor for jth MV, in engineering units.
$w_{i, j}^{u}$ — Tuning weight for jth MV at ith prediction horizon step (dimensionless).

The values n_u, p, $s_{j}^{u}$ , and $w_{i, j}^{u}$ are constant controller specifications. The controller receives u_j,target(k+i|k) values for the entire horizon. The controller uses the state observer to predict the plant outputs. Thus, J_u is a function of z_k only.

Manipulated Variable Move Suppression

Most applications prefer small MV adjustments (moves). An MPC constant uses the following scalar performance measure for manipulated variable move suppression:

$J_{Δ u} (z_{k}) = \sum_{j = 1}^{n_{u}} \sum_{i = 0}^{p - 1} {\frac{w_{i, j}^{Δ u}}{s_{j}^{u}} [u_{j} (k + i | k) - u_{j} (k + i - 1 | k)]}^{2} .$

Here,

k — Current control interval.
p — Prediction horizon (number of intervals).
n_u — Number of manipulated variables.
z_k — QP decision, given by:
$z_{k}^{T} = [\begin{matrix} u {(k | k)}^{T} & u {(k + 1 | k)}^{T} & \dots & u {(k + p - 1 | k)}^{T} & ε_{k} \end{matrix}] .$
$s_{j}^{u}$ — Scale factor for jth MV, in engineering units.
$w_{i, j}^{Δ u}$ — Tuning weight for jth MV movement at ith prediction horizon step (dimensionless).

The values n_u, p, $s_{j}^{u}$ , and $w_{i, j}^{Δ u}$ are constant controller specifications. u(k–1|k) = u(k–1), which are the known MVs from the previous control interval. J_Δu is a function of z_k only.

In addition, a control horizon m < p (or MV blocking) constrains certain MV moves to be zero.

Constraint Violation

In practice, constraint violations might be unavoidable. Soft constraints allow a feasible QP solution under such conditions. An MPC controller employs a dimensionless, nonnegative slack variable, ε_k, which quantifies the worst-case constraint violation. (See Constraints) The corresponding performance measure is:

$J_{ε} (z_{k}) = ρ_{ε} ε_{k}^{2} .$

Here,

z_k — QP decision, given by:
$z_{k}^{T} = [\begin{matrix} u {(k | k)}^{T} & u {(k + 1 | k)}^{T} & \dots & u {(k + p - 1 | k)}^{T} & ε_{k} \end{matrix}] .$
ε_k — Slack variable at control interval k (dimensionless).
ρ_ε — Constraint violation penalty weight (dimensionless).

Alternative Cost Function

You can elect to use the following alternative to the standard cost function:

$J (z_{k}) = \sum_{i = 0}^{p - 1} {[e_{y}^{T} (k + i) Q e_{y} (k + i)] + [e_{u}^{T} (k + i) R_{u} e_{u} (k + i)] + [Δ u^{T} (k + i) R_{Δ u} Δ u (k + i)]} + ρ_{ϵ} ε_{k}^{2} .$

Here, Q (n_y-by-n_y), R_u, and R_Δu (n_u-by-n_u) are positive-semi-definite weight matrices, and:

$\begin{matrix} e_{y} (i + k) = S_{y}^{- 1} [r (k + i + 1 | k) - y (k + i + 1 | k)] \\ e_{u} (i + k) = S_{u}^{- 1} [u_{t a r g e t} (k + i | k) - u (k + i | k)] \\ Δ u (k + i) = S_{u}^{- 1} [u (k + i | k) - u (k + i - 1 | k)] . \end{matrix}$

Also,

S_y — Diagonal matrix of plant output variable scale factors, in engineering units.
S_u — Diagonal matrix of MV scale factors in engineering units.
r(k+1|k) — n_y plant output reference values at the ith prediction horizon step, in engineering units.
y(k+1|k) — n_y plant outputs at the ith prediction horizon step, in engineering units.
z_k — QP decision, given by:
$z_{k}^{T} = [\begin{matrix} u {(k | k)}^{T} & u {(k + 1 | k)}^{T} & \dots & u {(k + p - 1 | k)}^{T} & ε_{k} \end{matrix}] .$
u_target(k+i|k) — n_u MV target values corresponding to u(k+i|k), in engineering units.

Output predictions use the state observer, as in the standard cost function.

The alternative cost function allows off-diagonal weighting, but requires the weights to be identical at each prediction horizon step.

The alternative and standard cost functions are identical if the following conditions hold:

The standard cost functions employs weights $w_{i, j}^{y}$ , $w_{i, j}^{u}$ , and $w_{i, j}^{Δ u}$ that are constant with respect to the index, i = 1:p.
The matrices Q, R_u, and R_Δu are diagonal with the squares of those weights as the diagonal elements.

Constraints

Certain constraints are implicit. For example, a control horizon m < p (or MV blocking) forces some MV increments to be zero, and the state observer used for plant output prediction is a set of implicit equality constraints. Explicit constraints that you can configure are described below.

Bounds on Plant Outputs, MVs, and MV Increments

The most common MPC constraints are bounds, as follows.

$\begin{matrix} \frac{y_{j, m i n} (i)}{s_{j}^{y}} - ε_{k} V_{j, m i n}^{y} (i) \leq \frac{y_{j} (k + i | k)}{s_{j}^{y}} \leq \frac{y_{j, m a x} (i)}{s_{j}^{y}} + ε_{k} V_{j, m a x}^{y} (i), i = 1 : p, j = 1 : n_{y} \\ \frac{u_{j, m i n} (i)}{s_{j}^{u}} - ε_{k} V_{j, m i n}^{u} (i) \leq \frac{u_{j} (k + i - 1 | k)}{s_{j}^{u}} \leq \frac{u_{j, m a x} (i)}{s_{j}^{u}} + ε_{k} V_{j, m a x}^{u} (i), i = 1 : p, j = 1 : n_{u} \\ \frac{Δ u_{j, m i n} (i)}{s_{j}^{u}} - ε_{k} V_{j, m i n}^{Δ u} (i) \leq \frac{Δ u_{j} (k + i - 1 | k)}{s_{j}^{u}} \leq \frac{Δ u_{j, m a x} (i)}{s_{j}^{u}} + ε_{k} V_{j, m a x}^{Δ u} (i), i = 1 : p, j = 1 : n_{u} . \end{matrix}$

Here, the V parameters (ECR values) are dimensionless controller constants analogous to the cost function weights but used for constraint softening (see Constraint Softening). Also,

ε_k — Scalar QP slack variable (dimensionless) used for constraint softening.
$s_{j}^{y}$ — Scale factor for jth plant output, in engineering units.
$s_{j}^{u}$ — Scale factor for jth MV, in engineering units.
y_j,min(i), y_j,max(i) — lower and upper bounds for jth plant output at ith prediction horizon step, in engineering units.
u_j,min(i), u_j,max(i) — lower and upper bounds for jth MV at ith prediction horizon step, in engineering units.
Δu_j,min(i), Δu_j,max(i) — lower and upper bounds for jth MV increment at ith prediction horizon step, in engineering units.

Except for the slack variable non-negativity condition, all of the above constraints are optional and are inactive by default (i.e., initialized with infinite limiting values). To include a bound constraint, you must specify a finite limit when you design the controller.

QP Matrices

This section describes the matrices associated with the model predictive control optimization problem described in Optimization Problem.

Prediction

Assume that the disturbance models described in Input Disturbance Model are unit gains; that is, d(k) = n_d(k) is white Gaussian noise. You can denote this problem as

$x \leftarrow [\begin{matrix} x \\ x_{d} \end{matrix}], A \leftarrow [\begin{matrix} A & B_{d} \bar{C} \\ 0 & \bar{A} \end{matrix}], B_{u} \leftarrow [\begin{matrix} B_{u} \\ 0 \end{matrix}], B_{v} \leftarrow [\begin{matrix} B v \\ 0 \end{matrix}], B_{d} \leftarrow [\begin{matrix} B_{d} \bar{D} \\ \bar{B} \end{matrix}], C \leftarrow [\begin{matrix} C & D_{d} \bar{C} \end{matrix}]$

Then, the prediction model is:

x(k+1) = Ax(k) +B_uu(k) +B_vv(k)+B_dn_d(k)

y(k) = Cx(k) +D_vv(k) +D_dn_d(k)

Next, consider the problem of predicting the future trajectories of the model performed at time k=0. Set n_d(i)=0 for all prediction instants i, and obtain

$y (i | 0) = C [A^{i} x (0) + \sum_{h = 0}^{i - 1} A^{i - 1} (B_{u} (u (- 1) + \sum_{j = 0}^{h} Δ u (j)) + B_{v} v (h))] + D_{v} v (i)$

This equation gives the solution

$[\begin{matrix} y (1) \\ \dots \\ y (p) \end{matrix}] = S_{x} x (0) + S_{u 1} u (- 1) + S_{u} [\begin{matrix} Δ u (0) \\ \dots \\ Δ u (p - 1) \end{matrix}] + H_{v} [\begin{matrix} v (0) \\ \dots \\ v (p) \end{matrix}]$

where

$\begin{array}{l} S_{x} = [\begin{matrix} C A \\ C A^{2} \\ \dots \\ C A^{p} \end{matrix}] \in ℜ^{p n_{y} \times n_{x}}, S_{u 1} = [\begin{matrix} C B_{u} \\ C B_{u} + C A B_{u} \\ \dots \\ \sum_{h = 0}^{p - 1} C A^{h} B_{u} \end{matrix}] \in ℜ^{p n_{y} \times n_{u}} \\ S_{u} = [\begin{matrix} C B_{u} & 0 & \dots & 0 \\ C B_{u} + C A B_{u} & C B_{u} & \dots & 0 \\ \dots & \dots & \dots & \dots \\ \sum_{h = 0}^{p - 1} C A^{h} B_{u} & \sum_{h = 0}^{p - 2} C A^{h} B_{u} & \dots & C B_{u} \end{matrix}] \in ℜ^{p n_{y} \times p n_{u}} \\ H_{v} = [\begin{matrix} C B_{v} & D_{v} & 0 & \dots & 0 \\ C A B_{v} & C B_{v} & D_{v} & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots \\ C A^{p - 1} B_{v} & C A^{p - 2} B_{v} & C A^{p - 3} B_{v} & \dots & D_{v} \end{matrix}] \in ℜ^{p n_{y} \times (p + 1) n_{v}} . \end{array}$

Optimization Variables

Let m be the number of free control moves, and let z= [z₀; ...; z_m–1]. Then,

$[\begin{matrix} Δ u (0) \\ \dots \\ Δ u (p - 1) \end{matrix}] = J_{M} [\begin{matrix} z_{0} \\ \dots \\ z_{m - 1} \end{matrix}]$

where J_M depends on the choice of blocking moves. Together with the slack variable ɛ, vectors z₀, ..., z_m–1 constitute the free optimization variables of the optimization problem. In the case of systems with a single manipulated variable, z₀, ..., z_m–1 are scalars.

Consider the blocking moves depicted in the following graph.

Blocking Moves: Inputs and Input Increments for moves = [2 3 2]

This graph corresponds to the choice moves=[2 3 2], or equivalently, u(0)=u(1), u(2)=u(3)=u(4), u(5)=u(6), Δ u(0)=z0, Δ u(2)=z1, Δ u(5)=z2, Δ u(1)=Δ u(3)=Δ u(4)=Δ u(6)=0.

Then, the corresponding matrix J_M is

$J_{M} = [\begin{matrix} I & 0 & 0 \\ 0 & 0 & 0 \\ 0 & I & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & I \\ 0 & 0 & 0 \end{matrix}]$

For more information on manipulated variable blocking, see Manipulated Variable Blocking.

Cost Function

Standard Form. The function to be optimized is

$\begin{array}{l} J (z, ε) = {([\begin{matrix} u (0) \\ \dots \\ u (p - 1) \end{matrix}] - [\begin{matrix} u_{t a r g e t} (0) \\ \dots \\ u_{t a r g e t} (p - 1) \end{matrix}])}^{T} W_{u}^{2} ([\begin{matrix} u (0) \\ \dots \\ u (p - 1) \end{matrix}] - [\begin{matrix} u_{t a r g e t} (0) \\ \dots \\ u_{t a r g e t} (p - 1) \end{matrix}]) + {[\begin{matrix} Δ u (0) \\ \dots \\ Δ u (p - 1) \end{matrix}]}^{T} W_{Δ u}^{2} [\begin{matrix} Δ u (0) \\ \dots \\ Δ u (p - 1) \end{matrix}] \\ + {([\begin{matrix} y (1) \\ \dots \\ y (p) \end{matrix}] - [\begin{matrix} r (1) \\ \dots \\ r (p) \end{matrix}])}^{T} W_{y}^{2} ([\begin{matrix} y (1) \\ \dots \\ y (p) \end{matrix}] - [\begin{matrix} r (1) \\ \dots \\ r (p) \end{matrix}]) + ρ_{ε} ε^{2} \end{array}$

where

\begin{matrix} W_{u} = diag (w_{0, 1}^{u}, w_{0, 2}^{u}, ..., w_{0, n_{u}}^{u}, ..., w_{p - 1, 1}^{u}, w_{p - 1, 2}^{u}, ..., w_{p - 1, n_{u}}^{u}) \\ W_{Δ u} = diag (w_{0, 1}^{Δ u}, w_{0, 2}^{Δ u}, ..., w_{0, n_{u}}^{Δ u}, ..., w_{p - 1, 1}^{Δ u}, w_{p - 1, 2}^{Δ u}, ..., w_{p - 1, n_{u}}^{Δ u}) \\ W_{y} = diag (w_{1, 1}^{y}, w_{1, 2}^{y}, ..., w_{1, n_{y}}^{y}, ..., w_{p, 1}^{y}, w_{p, 2}^{y}, ..., w_{p, n_{y}}^{y}) \end{matrix}

(1)

Finally, after substituting u(k), Δu(k), y(k), J(z) can be rewritten as

\begin{array}{l} J (z, ε) = ρ_{ε} ε^{2} + z^{T} K_{Δ u} z + 2 ({[\begin{matrix} r (1) \\ \dots \\ r (p) \end{matrix}]}^{T} K_{r} + {[\begin{matrix} v (0) \\ \dots \\ v (p) \end{matrix}]}^{T} K_{v} + u {(- 1)}^{T} K_{u} + {[\begin{matrix} u_{t a r g e t} (0) \\ \dots \\ u_{t a r g e t} (p - 1) \end{matrix}]}^{T} K_{u t} + x {(0)}^{T} K_{x}) z \\ + c_{y}^{T} W_{y} c_{y} + c_{u}^{T} W_{u} c_{u} \end{array}

(2)

where

$\begin{array}{l} c_{y} = S_{x} x (0) + S_{u 1} u (- 1) + H_{v} [\begin{matrix} v (0) \\ \dots \\ v (p) \end{matrix}] - [\begin{matrix} r (1) \\ \dots \\ r (p) \end{matrix}] \\ c_{u} = [\begin{matrix} I_{1} \\ \dots \\ I_{p} \end{matrix}] u (- 1) - [\begin{matrix} u_{t a r g e t} (0) \\ \dots \\ u_{t a r g e t} (p - 1) \end{matrix}] \end{array}$

Here, I₁ = … = I_p are identity matrices of size n_u.

Note

You may want the QP problem to remain strictly convex. If the condition number of the Hessian matrix K_Δ_U is larger than 10¹², add the quantity 10*sqrt(eps) on each diagonal term. You can use this solution only when all input rates are unpenalized (W^Δ^u=0) (see the Weights property of the mpc object).

Alternative Cost Function. If you are using the alternative cost function shown in Alternative Cost Function, then Equation 1 is replaced by the following:

\begin{array}{l} W_{u} = blkdiag (R_{u}, ..., R_{u}) \\ W_{Δ u} = blkdiag (R_{Δ u}, ..., R_{Δ u}) \\ W_{y} = blkdiag (Q, ..., Q) \end{array}

(3)

In this case, the block-diagonal matrices repeat p times, for example, once for each step in the prediction horizon.

You also have the option to use a combination of the standard and alternative forms. For more information, see the Weights property of the mpc object.

Constraints

Next, consider the limits on inputs, input increments, and outputs along with the constraint ɛ≥ 0.

$[\begin{matrix} y_{\min} (1) - ε V_{\min}^{y} (1) \\ \dots \\ y_{\min} (p) - ε V_{\min}^{y} (p) \\ u_{\min} (0) - ε V_{\min}^{u} (0) \\ \dots \\ u_{\min} (p - 1) - ε V_{\min}^{u} (p - 1) \\ Δ u_{\min} (0) - ε V_{\min}^{Δ u} (0) \\ \dots \\ Δ u_{\min} (p - 1) - ε V_{\min}^{Δ u} (p - 1) \end{matrix}] \leq [\begin{matrix} y (1) \\ \dots \\ y (p) \\ u (0) \\ \dots \\ u (p - 1) \\ Δ u (0) \\ \dots \\ Δ u (p - 1) \end{matrix}] \leq [\begin{matrix} y_{\max} (1) + ε V_{\max}^{y} (1) \\ \dots \\ y_{\max} (p) + ε V_{\max}^{y} (p) \\ u_{\max} (0) + ε V_{\max}^{u} (0) \\ \dots \\ u_{\max} (p - 1) + ε V_{\max}^{u} (p - 1) \\ Δ u_{\max} (0) + ε V_{\max}^{Δ u} (0) \\ \dots \\ Δ u_{\max} (p - 1) + ε V_{\max}^{Δ u} (p - 1) \end{matrix}]$

Note

To reduce computational effort, the controller automatically eliminates extraneous constraints, such as infinite bounds. Thus, the constraint set used in real time may be much smaller than that suggested in this section.

Similar to what you did for the cost function, you can substitute u(k), Δu(k), y(k), and obtain

M_{z} z + M_{ε} ε \leq M_{\lim} + M_{v} [\begin{matrix} v (0) \\ \dots \\ v (p) \end{matrix}] + M_{u} u (- 1) + M_{x} x (0)

(4)

In this case, matrices M_z, M_ɛ, M_lim, M_v, M_u, and M_x are obtained from the upper and lower bounds and ECR values.

Unconstrained Model Predictive Control

The optimal solution is computed analytically

$z * = - K_{Δ u}^{- 1} {({[\begin{matrix} r (1) \\ \dots \\ r (p) \end{matrix}]}^{T} K_{r} + [\begin{matrix} v (0) \\ \dots \\ v (p) \end{matrix}] K_{v} + u {(- 1)}^{T} K_{u} + {[\begin{matrix} u_{t a r g e t} (0) \\ \dots \\ u_{t a r g e t} (p - 1) \end{matrix}]}^{T} K_{u t} + x {(0)}^{T} K_{x})}^{T}$

and the model predictive controller sets Δu(k)=z*₀, u(k)=u(k–1)+Δu(k).

Documentation