Model Predictive Control

Question

I have a few confusions about Model Predictive Control (MPC). Since they are all minor questions related to the same category, I ask them under one topic.

In an article, the cost function is defined as: $$J(t)=\sum_{j=1}^{N_p}\delta(j) ( y(t+j|t) -ref(t+j) )^2 + \sum_{j=0}^{N_c-1}\lambda(j) u(t+j)^2+ \sum_{j=0}^{N_c-1}\gamma(j) \Delta u(t+j)^2$$

where, y is output, ref is the reference and u is the modified input, $N_C$ control horizon and $N_P$ the prediction window size.

1- First of all, in MPC, how to know reference in future? Even though we know the model, we do not know the future input. Is input considered to be constant? Unfortunately, I found no article explaining that.

2- In MPC, are $\delta(j), \lambda(j), \gamma(j)$ uniform? if they reduce over $j$, How do they usually reduce? exponentially?

3- What is the algorithm to guess $u(t+j)$ ? Are they guessed randomly? or u(t+j) depends on u(t+j-1) while initial guess (I considered their selection is done through Genetic Algorithm but I do not think it is the right algorithm for MPC).

4- Why $y(t+j|t)$ is predicted in continuous form while the manipulated input is discrete? Wouldn't it be better to be continuous too? At least by connection of dots with a simple line? (see the figure below)

5- How is the above equation simplified to: $$ J=\frac{1}2 \Delta U^T H \Delta U +\Delta U^T F $$ What happened to $y$? is $y$ a part of new matrix $U$? (article: Wang, Liuping. Model predictive control system design and implementation using MATLAB®. springer, 2009.)

.

Answer 1

1. In some cases you know the reference in advance (a robot supposed to follow a predefined trajectory). In some cases you don't (a pilot giving reference inputs) and you have to make a reasonable assumption (e.g., reference stays the same as current value)

2. Most often constant. If you want to use something else, you are free to come up with what ever you want (expontial decay, geometric decay etc etc)

3. Standard MPC boils down to simple convex quadratic programming (QP) which is a well established field with many solvers.

4. What is continuous? MPC isalmost always defined using the discrete-time model $x_{k+1} = Ax_k + Bu_k, y_k = Cx_k + D u_k$

5. By using the dicsrete-time model above and writing everything using vectors and matrices. $y$ can be eliminated if one wants as it is linear function of $x$ which is a linear function of the initial state and the control sequence. You will find this in any basic book on MPC.

Answer 2

In addition to the answer of @Johan at point 5), note I use a different cost function, with which you are probably more familiar with

Define the cost function $\mathbf{V}(\mathbf{x},\mathbf{U})$ as

$$ \mathbf{V}(\mathbf{x},\mathbf{U}) := \mathbf{x}_N^\mathrm{T} \mathbf{P} \mathbf{x}_N + \sum_{i=0}^{N-1}\left( \mathbf{x}_i^\mathrm{T} \mathbf{Q} \mathbf{x}_i + \mathbf{u}_i^\mathrm{T} \mathbf{R} \mathbf{u}_i\right) $$

Herein $N$ is the control horizon, $\mathbf{P}$ is the terminal state weight, $\mathbf{Q}$ is the state weight and $\mathbf{R}$ is the input weight. Remark: $\mathbf{P}$, $\mathbf{Q}$ and $\mathbf{R}$ are positive semidefinite. Furthermore, $\mathbf{U}$ is the stacked input vector, that is,

$$ \mathbf{U} := \begin{bmatrix} \mathbf{u}_0^T & \mathbf{u}_1^T & \ldots & \mathbf{u}_{N-1}^T \end{bmatrix}^T $$

The optimal control $\mathbf{U}^*$ input is given by

$$ \mathbf{U}^* = \arg \mathbf{V}^*(\mathbf{x},\mathbf{U}) = \arg \min_{\mathbf{U}} \mathbf{V}(\mathbf{x},\mathbf{U}) $$

The sum in the cost function can be rewritten in a matrix equation

$$ \begin{gathered} \mathbf{V}(\mathbf{x},\mathbf{U}) = \mathbf{x}_0^\mathrm{T} \mathbf{Q} \mathbf{x}_0 + \begin{bmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_N \end{bmatrix}^\mathrm{T} \begin{bmatrix} \mathbf{Q} & 0 & \cdots & 0 \\ 0 & \mathbf{Q} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathbf{P} \end{bmatrix} \begin{bmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_N \end{bmatrix} + \begin{bmatrix} \mathbf{u}_0 \\ \mathbf{u}_1 \\ \vdots \\ \mathbf{u}_{N-1} \end{bmatrix}^\mathrm{T} \begin{bmatrix} \mathbf{R} & 0 & \cdots & 0 \\ 0 & \mathbf{R} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathbf{R} \end{bmatrix} \begin{bmatrix} \mathbf{u}_0 \\ \mathbf{u}_1 \\ \vdots \\ \mathbf{u}_{N-1} \end{bmatrix} \\ \Leftrightarrow \\ \mathbf{V}(\mathbf{x},\mathbf{U}) = \mathbf{x}^\mathrm{T}\mathbf{Q}\mathbf{x} + \mathbf{X}^\mathrm{T}\mathbf{\Omega}\mathbf{X} + \mathbf{U}^\mathrm{T}\mathbf{\Psi}\mathbf{U} \end{gathered} $$

Remark: $\mathbf{x} = \mathbf{x}_0$. Furthermore, $\mathbf{P} \succeq 0 \wedge \mathbf{Q} \succeq 0$ implies that $\mathbf{\Omega} \succeq 0$ and $\mathbf{R} \succ 0$ implies that $\mathbf{\Psi} \succ 0$. The cost function now still depends on the future state $\mathbf{X}$. To eliminate $\mathbf{X}$ the prediction matrices $\mathbf{\Phi}$ and $\mathbf{\Gamma}$ are computed such that the stacked state vector is a function of the initial state and the stacked input vector i.e. $\mathbf{X} = \mathbf{\Phi}\mathbf{x} + \mathbf{\Gamma}\mathbf{U}$.

$$ \mathbf{\Phi} = \begin{bmatrix} \mathbf{A} \\ \mathbf{A}^2 \\ \vdots \\ \mathbf{A}^N \end{bmatrix} \qquad \mathbf{\Gamma} = \begin{bmatrix} \mathbf{B} & 0 & \cdots & 0 \\ \mathbf{A}\mathbf{B} & \mathbf{B} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{A}^{N-1}\mathbf{B} & \mathbf{A}^{N-2}\mathbf{B} & \cdots & \mathbf{B} \end{bmatrix} $$

The cost function $\mathbf{V}(\mathbf{x},\mathbf{U})$ can now be written as a function of the initial state $\mathbf{x}$ and the stacked input vector $\mathbf{U}$

$$ \mathbf{V}(\mathbf{x},\mathbf{U}) = \frac{1}{2} \mathbf{U}^\mathrm{T} \mathbf{G} \mathbf{U} + \mathbf{U}^\mathrm{T} \mathbf{F}\mathbf{x} + \mathbf{x}^\mathrm{T} \left( \mathbf{Q} + \mathbf{\Phi}^\mathrm{T} \mathbf{\Omega} \mathbf{\Phi} \right) \mathbf{x} $$

Where $\mathbf{G} = 2\left(\mathbf{\Psi} + \mathbf{\Gamma}^\mathrm{T} \mathbf{\Omega} \mathbf{\Gamma} \right)$, $\mathbf{F} = 2 \mathbf{\Gamma}^\mathrm{T} \mathbf{\Omega} \mathbf{\Phi}$ and $\mathbf{G} \succ 0$. The last term of the newly derived cost function is independent on $\mathbf{U}$. Therefore, it can be neglected when computing the optimal input sequence $\mathbf{U}^*(\mathbf{x})$

Now lets add some input and state constraints

$$ \mathbf{u}_\mathrm{min} \leq \mathbf{u}_k \leq \mathbf{u}_\mathrm{max} \wedge \mathbf{x}_\mathrm{min} \leq \mathbf{x}_k \leq \mathbf{x}_\mathrm{max} $$

This can be rewritten to

$$ \begin{bmatrix} 0 \\ 0 \\ -\mathbf{I}_n \\ \mathbf{I}_n \end{bmatrix} \mathbf{x}_k + \begin{bmatrix} -\mathbf{I}_m \\ \mathbf{I}_m \\ 0 \\ 0 \end{bmatrix} \mathbf{u}_k \leq \begin{bmatrix} -\mathbf{u}_\mathrm{min} \\ \mathbf{u}_\mathrm{max} \\ -\mathbf{x}_\mathrm{min} \\ \mathbf{x}_\mathrm{max} \end{bmatrix} \Leftrightarrow \mathbf{M}_k\mathbf{x}_k + \mathbf{E}_k\mathbf{u}_k \leq \mathbf{b}_k $$

Now note that this is only valid for $k = 0,1 \ldots, N-1$ so not up to the control horizon $N$. As such we add a terminal constraint for $N$ on the state.

$$ \begin{bmatrix} -\mathbf{I}_{n} \\ \mathbf{I}_{n} \\ \end{bmatrix} \mathbf{x}_N \leq \begin{bmatrix} -\mathbf{x}_\mathrm{min} \\ \mathbf{x}_\mathrm{max} \\ \end{bmatrix} \Rightarrow \mathbf{M}_N\mathbf{x}_N \leq \mathbf{b}_N $$

All constraints can now be combined to

$$ \begin{gathered} \begin{bmatrix} \mathbf{M}_0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \mathbf{x}_0 + \begin{bmatrix} 0 & \cdots & 0 \\ \mathbf{M}_1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \mathbf{M}_N \end{bmatrix} \begin{bmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_N \end{bmatrix} + \begin{bmatrix} \mathbf{E}_0 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \mathbf{E}_{N-1} \\ 0 & \cdots & 0 \end{bmatrix} \begin{bmatrix} \mathbf{u}_0 \\ \mathbf{u}_1 \\ \vdots \\ \mathbf{u}_N \end{bmatrix} \leq \begin{bmatrix} \mathbf{b}_0 \\ \mathbf{b}_1 \\ \vdots \\ \mathbf{b}_N \end{bmatrix} \\ \Leftrightarrow \\ \mathbf{\mathcal{D}}\mathbf{x} + \mathbf{\mathcal{M}}\mathbf{X} + \mathbf{\mathcal{E}}\mathbf{U} \leq \mathbf{c} \end{gathered} $$

Again this equation still depends on $\mathbf{X}$, which we eliminate using the same procedure as before. As a result, the above equation can be written as

$$\mathbf{J}\mathbf{U} \leq \mathbf{c} + \mathbf{W}\mathbf{x}$$

With $\mathbf{J} = \mathbf{\mathcal{M}}\mathbf{\Gamma} + \mathbf{\mathcal{E}}$ and $\mathbf{W} = -\mathbf{\mathcal{D}} - \mathbf{\mathcal{M}}\mathbf{\Phi}$. Now the constraints only depend on the stacked input vector and the initial state. The constrained MPC problem can then be formulated as

$$ \arg \underset{\mathbf{U}}{\min}\; \mathbf{V}(\mathbf{x},\mathbf{U}) \quad \text{s.t.} \quad \mathbf{J}\mathbf{U} \leq \mathbf{c} + \mathbf{W}\mathbf{x} $$

Now if you use a QP solver you will have the following QP problem

$$ \begin{aligned} \mathbf{U}^*(\mathbf{x}) = \arg \underset{\mathbf{U}}{\min}\;\;& \frac{1}{2}\mathbf{U}^\mathrm{T} \mathbf{G} \mathbf{U} + \mathbf{U}^\mathrm{T} \mathbf{F} \mathbf{x} \\[1ex] \text{subject to:}\;\;&\mathbf{x}_0 = \mathbf{x} \hspace{3cm} \text{Measurement} \nonumber \\ &\mathbf{J}\mathbf{U} \leq \mathbf{c} + \mathbf{W}\mathbf{x} \hspace{1.82cm} \text{Constraints} \nonumber \\ &\mathbf{Q} \succ 0, \mathbf{R} \succ 0 \hspace{1.99cm} \text{Cost weights} \end{aligned} $$

Now if you want to solve these kind of typical MPC problems, I would strongly advice to use the MPT toolbox, http://people.ee.ethz.ch/~mpt/3/, an example http://control.ee.ethz.ch/~mpt/3/UI/RegulationProblem. Furthermore, the slides from this workshop may help you as well https://www.kirp.chtf.stuba.sk/pc11/data/workshops/mpc/