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 variables vector, 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 the jth plant output at the ith prediction horizon step, in engineering units
r_j(k+i|k) — Reference value for the jth plant output at the ith prediction horizon step, in engineering units
$s_{j}^{y}$ — Scale factor for the jth plant output, in engineering units.
$w_{i, j}^{y}$ — Tuning weight for the jth plant output at the ith prediction horizon step (dimensionless).
ε_k — Slack variable at control interval k (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.

The problem formulation used here, in which the decision variables of the optimization problem are the are just the manipulated variables (plus a slack variable), is called a dense formulation. it is also known as a sequential or single-shooting formulation.

A different formulation is possible in which the decision variables of the optimization problem also include all the intermediate state variables. In this case, equality constraints of the type x(k+1) = Ax(k) + Bu(k) must be enforced. This formulation results in a larger problem with many equality constraints, but a more regular sparse structure, an easier-to-compute cost function, and better numerical conditioning for large horizons. This alternative formulation is referred to as sparse, (but also simultaneous or multiple-shooting) approach.

The default problem formulation used in Model Predictive Control Toolbox™ for linear MPC problems is the dense formulation, because it can have a smaller memory footprint (if a generic QP solver is used) and tends to be more efficient for stable linear problems in which the prediction horizon is small and when a considerable amount of Manipulated Variable Blocking is used.

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 variables vector, 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 the jth MV at the ith prediction horizon step, in engineering units
$s_{j}^{u}$ — Scale factor for the jth MV, in engineering units
$w_{i, j}^{u}$ — Tuning weight for the jth MV at the ith prediction horizon step (dimensionless)
ε_k — Slack variable at control interval k (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 variables vector, 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 the jth MV, in engineering units
$w_{i, j}^{Δ u}$ — Tuning weight for the jth MV movement at the ith prediction horizon step (dimensionless)
ε_k — Slack variable at control interval k (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 variables vector, 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 the control interval k (dimensionless).
ρ_ε — Constraint violation penalty weight (dimensionless). This weight corresponds to the Weights.ECR property of an mpc object.

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) = L_{S y}^{- 1} [r (k + i + 1 | k) - y (k + i + 1 | k)] \\ e_{u} (i + k) = L_{S u}^{- 1} [u_{t a r g e t} (k + i | k) - u (k + i | k)] \\ Δ u (k + i) = L_{S u}^{- 1} [u (k + i | k) - u (k + i - 1 | k)] \end{matrix}$

Also,

L_Sy — Diagonal matrix of plant output variable scale factors s_y, in engineering units
L_Su — Diagonal matrix of MV scale factors s_u, 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 variables vector, 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, and correspond to the MinECR, MaxECR, RateMinECR and RateMaxECR vectors stored in the OutputVariables and ManipulatedVariables properties of an mpc object. For more information, see Specify Constraints.

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.

Mixed Input/Output Linear Constraints

You can also introduce mixed input/output linear constraints.

The general form of the mixed input/output constraints is:

Eu(k + j) + Fy(k + j) + Sv(k + j) ≤ G + εV

Here, j = 0,...,p, and:

p is the prediction horizon.
k is the current time index.
u is a column vector manipulated variables.
y is a column vector of all plant output variables.
v is a column vector of measured disturbance variables.
ε is a scalar slack variable used for constraint softening (as in Standard Cost Function).
E, F, G, V, and S are constant matrices.

For more information, see Constraints on Linear Combinations of Inputs and Outputs, setconstraint, getconstraint and Update Mixed Input/Output Constraints at Run Time.

QP Matrices

This section describes the matrices associated with the model predictive control optimization problem described in Optimization Problem. Note that the vector of decision variables used in this section (and in the actual implementation of the software) is:

$z_{k}^{T} = [\begin{matrix} Δ u {(k | k)}^{T} & Δ u {(k + 1 | k)}^{T} & \dots & Δ u {(k + p - 1 | k)}^{T} & ε_{k} \end{matrix}] .$

where

$Δ u (k + i) = S_{u}^{- 1} [u (k + i | k) - u (k + i - 1 | k)] .$

This vector of decision variables is mathematically equivalent to the one described in the previous sections, but, from an implementation perspective, has advantages for formulations that allow for weights and constraints on increments of manipulated variables.

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} = L_{S u}^{- 1} 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} = L_{S u}^{- 1} 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} = L_{S y}^{- 1} 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)

and L_Sy and L_Su are the diagonal matrices or outputs and MV scaling factors, respectively, in engineering units.

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} = L_{S u}^{- 1} blkdiag (R_{u}, ..., R_{u}) \\ W_{Δ u} = L_{S u}^{- 1} blkdiag (R_{Δ u}, ..., R_{Δ u}) \\ W_{y} = L_{S y}^{- 1} 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).