Control theory: how do you initialize input for a model predictive controller?

Although I think this is a pretty popular algorithm, I posted this question on electronics stacksexchange but received no response. I am posting it here to see if there are many signal/control experts who can answer my very simple question

On these slides, there is a description of MPC algorithm. I am confused as to how the input $u$ is initialized and how its values are known in the future.

1. How is $u$ initialized at the beginning?

In the very beginning, we need all the future values of $u$ in order to predict the value of $x$. I am fully aware of how $x$ is predicted $N$ steps into the future.

• But where do you obtain $u_k, u_{k+1}, \ldots, u_{k+N-1}$ to predict the states in the first place?
• For instance how do you know $u_{k+1}$?
2. How to get the next $u$ values after receding horizon step?

After obtaining the optimal $u^*$ from the optimization program, we keep the first element of $u^*$ and discards the rest. This allows us to produce our output $y$.

• But how are the subsequent $u$ generated in order to feed the next iteration of the program?
• Again, we need every single value $u_k, u_{k+1}, \ldots, u_{k+N-1}$, but for instance how do you know $u_{k+1}$? (This is the same question as question 1 I suppose)

All in all, I am not sure how we know the future of our input $u$ in order to predict $x$. Can someone please enlighten me as to how this is done in simulation or in practice?

This algorithm in general tries to solve an optimization problem each step, defined as,

\begin{aligned} & \underset{\textbf{u}}{\text{minimize}} & & \sum_{i\ =\ 0}^{N-1}\left[x^T\!(k+1+i)\, Q\, x(k+1+i) + u^T\!(k+i)\, R\, u(k+i)\right] \\ & \text{subject to} & & x(k+1+i) = f(x(k+i),\ u(k+i)) \qquad \forall\ i = 0, \ldots, N \\ & & & g(x(k+i),\ u(k+i)) \leq 0 \qquad \qquad \qquad \ \ \ \forall\ i = 0, \ldots, N \end{aligned} \tag{1}

where,

$$\textbf{u} = \begin{bmatrix} u(k) & u(k+1) & \cdots & u(k+N-1) \end{bmatrix}^T. \tag{2}$$

This means that the algorithm tries to find the values for $\textbf{u}$ which minimize that quadratic summation expression, while still satisfying the update equality constraints and the inequality constraints (which returns a vector of length $m$). Here each $x(k)$ is the error between between the actual/predicted state and the desired reference state at $k$.

Depending on the functions $f(x,u)$ and $g(x,u)$ this problem can be hard to solve or even to prove that you found the global optimum and often require an initial guess. For example a generic solver for these kind of problem in Matlab can be solved with fmincon.

However if $f(x,u)$ and $g(x,u)$ are linear functions, then this problem becomes "a lot" easier to solve. A side note: if $g(x,u)$ is not linear, but is convex, then you can approximate it with linear inequalities. In Matlab the function quadprog can be used, since the const function is already in a quadratic form.

In chapter 2 of the paper it is shown how you can simplify the problem when this is the case. Namely all equality constraints and future states on the system can be removed. It also introduces $\bar Q$, which allows you to calculate the cost function from $N$ to infinity, assuming that you are no long constraint by the inequality constraints. Such an infinite horizon in the cost function helps ensuring stability of the system. Finding this $\bar Q$ involves solving the Riccati equation.

For example for a problem with no inequality constraints and only linear equality constraints, such that,

$$x(k+1+i) = A\, x(k+i) + B\, u(k+i) \qquad \forall\ i = 0, \ldots, N , \tag{3}$$

with $x(k) = x_0$, you can find all other $N$ states of the system as a function of $\textbf{u}$ with,

$$\begin{bmatrix} x(k+1) \\ x(k+2) \\ \vdots \\ x(k+N) \end{bmatrix} = \underbrace{\begin{bmatrix} A \\ A^2 \\ \vdots \\ A^N \end{bmatrix}}_{\textbf{A}} x_0 + \underbrace{\begin{bmatrix} B & 0 & \cdots & 0 \\ AB & B & \cdots & 0 \\ \vdots & \vdots & \ddots \\ A^{N-1}B & A^{N-2}B & \cdots & B \end{bmatrix}}_{\textbf{B}} \textbf{u}. \tag{4}$$

The cost function from equation $(1)$, denoted as $J(k)$ in your linked paper, can then also be written as,

$$J(k) = \left[\textbf{A} x_0\! + \textbf{B} \textbf{u}\right]^T\! \underbrace{\begin{bmatrix} Q & 0 & \cdots & 0 \\ 0 & Q & & \vdots \\ \vdots & & \ddots & 0 \\ 0 & \cdots & 0 & Q \end{bmatrix}}_{\textbf{Q}} \left[\textbf{A} x_0\! + \textbf{B} \textbf{u}\right] + \textbf{u}^T\! \underbrace{\begin{bmatrix} R & 0 & \cdots & 0 \\ 0 & R & & \vdots \\ \vdots & & \ddots & 0 \\ 0 & \cdots & 0 & R \end{bmatrix}}_{\textbf{R}} \textbf{u}. \tag{5}$$

After expanding the brackets, then equation $(5)$ can also be written as,

$$J(k) = x_0^T\textbf{A}^T\textbf{Q}\textbf{A}x_0 + x_0^T\textbf{A}^T\textbf{Q}\textbf{B}\textbf{u} + \textbf{u}^T\textbf{B}^T\textbf{Q}\textbf{A}x_0 + \textbf{u}^T \left[\textbf{B}^T\textbf{Q}\textbf{B} + \textbf{R}\right] \textbf{u}. \tag{6}$$

The first term of equation $(6)$ does not depend on $\textbf{u}$ and thus will be a constant. The second term is the transpose of the third term because $\textbf{Q}$ is symmetric, since their results are also a constant (a one by one matrix) they are equal and can be added together. Equation $(6)$ can therefore be simplified even further to,

$$J(k) = \underbrace{x_0^T\textbf{A}^T\textbf{Q}\textbf{A}x_0}_c + 2\,\underbrace{x_0^T\textbf{A}^T\textbf{Q}\textbf{B}}_{f}\,\textbf{u} + \textbf{u}^T\,\underbrace{\left[\textbf{B}^T\textbf{Q}\textbf{B} + \textbf{R}\right]}_{H}\,\textbf{u}. \tag{7}$$

Because it is assumed that $Q$ positive semi-definite and $R$ is positive definite, therefore the minimum of $J(k)$ can be found by taking its derivative with respect to $\textbf{u}$ and equating it to zero. This allows you to get the following solution for $\textbf{u}$, similar to equation $(2.6)$ of the paper,

$$\textbf{u} = -H^{-1}f. \tag{8}$$

So using the model from equation $(3)$ allows you to express the cost function only in terms of the future inputs. The solutions of problems with inequality constraints can be found in a similar way, but may require numerical algorithms however that is an entire field of study by itself.

• Hm I am still not sure because to predict the future values $x(k+i)$, for a linear system of the type $x(k+1) = Ax(k)+Bu(k)$, you find $x(k+i)$ using $x(k+i) = A^ix(k) + [A^{i-1}B | \ldots | B] *[u(k), u(k+1), \ldots, u(k+i-1)]^T$. (I probably have an index or two off) But you have no idea what the values of $u(k+1), \ldots$ are. The program doesn't solve for $x$, and instead takes it as an input. How can you obtain the full $x$ vector for the optimization without knowing the future values of $u$ before hand? – Carlos - the Mongoose - Danger Aug 24 '16 at 2:46
• I am confused because the first tutorial notes I found said you must have the full vector $x$ before the optimization step...See page 4 of eng.ox.ac.uk/~conmrc/mpc/mpc_lec1.pdf. It doesn't explain where $u$ comes from and it is written in a way that implies you need to have $u$ before optimizing – Carlos - the Mongoose - Danger Aug 24 '16 at 2:48
• @BeachedWhale I added an explanation for how to solve the unconstrained optimization problem, so I hope this helps you understand the underlying mechanisms for finding the future inputs (and indirectly also the future sates of the system). – fibonatic Aug 25 '16 at 17:07