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.
