plain##1 ##2 ##1 ##2 (##3)

# Numerical Discrete-Time Implementation
of Continuous-Time Linear-Quadratic Model Predictive Control

Zhanhao Zhang, Anders Hilmar Damm Christensen, Steen Hørsholt, John Bagterp
Jørgensen

###### Abstract

This study presents the design, discretization and implementation of the
continuous-time linear-quadratic model predictive control (CT-LMPC). The
control model of the CT-LMPC is parameterized as transfer functions with time
delays, and they are separated into deterministic and stochastic parts for
relevant control and filtering algorithms. We formulate time-delay, finite-
horizon CT linear-quadratic optimal control problems (LQ-OCPs) for the CT-
LMPC. By assuming piece-wise constant inputs and constraints, we present the
numerical discretization of the proposed LQ-OCPs and show how to convert the
discrete-time (DT) equivalent into a standard quadratic program. The
performance of the CT-LMPC is compared with the conventional DT-LMPC
algorithm. Our numerical experiments show that, under fixed tunning
parameters, the CT-LMPC shows better closed-loop performance as the sampling
time increases than the conventional DT-LMPC.

##  I Introduction

Consider the following continuous-time (CT) linear-quadratic optimal control
problem (LQ-OCP) for the linear model predictive control (LMPC)

$\displaystyle\min\>\>$
$\displaystyle\phi=\int_{t_{0}}^{t_{0}+T}l(\tilde{z}(t))dt$ (1a)
$\displaystyle s.t.$ $\displaystyle x(t_{0})=\hat{x}_{0},$ (1b) $\displaystyle
u(t)=u_{k},\quad$ $\displaystyle t_{k}\leq t<t_{k+1},\,k\in\mathcal{N},$ (1c)
$\displaystyle\dot{x}(t)=A_{c}x(t)+B_{c}u(t),\quad$ $\displaystyle t_{0}\leq
t<t_{0}+T,$ (1d) $\displaystyle z(t)=C_{c}x(t)+D_{c}u(t),\,$ $\displaystyle
t_{0}\leq t<t_{0}+T,$ (1e) $\displaystyle\bar{z}(t)=\bar{z}_{k},$
$\displaystyle t_{k}\leq t<t_{k+1},\,k\in\mathcal{N},$ (1f)
$\displaystyle\tilde{z}(t)=z(t)-\bar{z}(t),$ $\displaystyle t_{0}\leq
t<t_{0}+T,$ (1g)

with the stage cost function

$\begin{split}l_{c}(\tilde{z}(t))&=\frac{1}{2}\left\|W_{z}\tilde{z}(t)\right\|_{2}^{2}=\frac{1}{2}\left(\tilde{z}(t)^{\prime}Q_{c}\tilde{z}(t)\right),\end{split}$
(2)

where $x\in\mathbb{R}^{n_{x}\times 1}$ is the state,
$u\in\mathbb{R}^{n_{u}\times 1}$ is the input, $z\in\mathbb{R}^{n_{z}\times
1}$ and $\bar{z}\in\mathbb{R}^{n_{z}\times 1}$ are the output and its
reference. $Q_{c}=W_{z}^{\prime}W_{z}$ is semi-positive definite weight
matrix. $T=NT_{s}$ is the control horizon and $\mathcal{N}=0,1,\ldots,N-1$.

In (1), we are assuming piece-wise constant inputs (1c) and reference (1f)
since most control systems in real-world scenarios operate on digital
platforms. The corresponding discrete-time (DT) equivalent is

$\displaystyle\min_{x,u}\quad$
$\displaystyle\phi=\sum_{k\in\mathcal{N}}l_{k}(x_{k},u_{k})$ (3a)
$\displaystyle s.t.$ $\displaystyle x_{0}=\hat{x}_{0},$ (3b) $\displaystyle
x_{k+1}=Ax_{k}+Bu_{k},\quad$ $\displaystyle k\in\mathcal{N},$ (3c)

with the stage costs

$l_{k}(x_{k},u_{k})=\frac{1}{2}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}^{\prime}Q\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}+q_{k}^{\prime}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}+\rho_{k},\quad k\in\mathcal{N},$ (4)

where $Q$ is a symmetric positive semi-definite matrix and

$Q=\begin{bmatrix}Q_{xx}&Q_{xu}\\\ Q_{ux}&Q_{uu}\end{bmatrix},\qquad
q_{k}=\begin{bmatrix}q_{x,k}\\\ q_{u,k}\end{bmatrix}.$ (5)

The LQ-OCP is considered one of the fundamental problems in optimal control
theory. This class of problems capture plenty of attention due to their
simplicity and mathematical tractability. There is abundant research on the
solution and discretization methods of optimal control problems described in
(1) and (3). Most studies consider piece-wise constant inputs and use
numerical discretization methods (e.g., Euler and Runge-Kutta methods) to
obtain the DT equivalent and derive the solution [1, 9, 4, 6, 20]. In [2, 3,
5], authors explored the discrete equivalent of the LQ-OCP (1) and described
the analytic expressions of the desired equivalent, discrete weighting
matrices $Q$ and $q_{k}$ described in (5).

On the other hand, model predictive control, which might be the most
successful advanced control algorithm in the industry, has a strong connection
with the LQ-OCP. This relationship arises from using the LQ-OCP as the
underlying optimization problem within the MPC framework. However, due to the
consideration of computational efficiency and ease of implementation, almost
all LMPC applications implemented in practice use DT models and diagonal
weight matrices rather than the discrete weight matrix $Q$ described in (5).
In principle, non-diagonal weight matrices can be used. However, this is
almost infeasible in practice as it will significantly increase the complexity
of tuning. The previous studies [17, 16] proposed a novel computational
procedure for input-constrained CT linear-quadratic regulator problem (LQR)
with an infinite horizon. They pointed out that the discrete weight matrix $Q$
can be computed with one exponential matrix instead of standard numerical
discretization methods, and the numerical experiments showed that the CT LQR
offers significant advantages over the DT MPC. Some other interesting results
about the CT LQR or CT MPC can be found in [12, 21, 13].

However, the existing literature has not investigated the case with time
delays, which is common in real-world scenarios (e.g., the system model is a
MIMO transfer function with time delays). Therefore, in this paper, we would
like to investigate the CT linear-quadratic model predictive control (CT-LMPC)
and compare it with the standard DT-LMPC. This paper focuses on the following
problems:

  * 1.

The numerical discretization of CT LQ-OCPs subject to linear time-delay
systems

  * 2.

The numerical discrete-time implementation of CT-LMPC

This paper is organized as follows. Section II introduces the control model
and different objective functions of the CT-LMPC and present their
discretization and implementation. Section III presents numerical experiments
comparing and testing the proposed CT-LMPC and conventional DT-LMPC. The
conclusions are given in Section IV.

##  II Discretization of CT-LMPC

This section describes the discretization of linear time-delay models and the
LQ-OCPs for the CT-LMPC.

###  II-A Discretization of linear time-delay models

Consider a MIMO input-output model

$\displaystyle\boldsymbol{Z}(s)$
$\displaystyle=G(s)U(s)+H(s)\boldsymbol{W}(s),$ (6a)
$\displaystyle\boldsymbol{Y}(s)$
$\displaystyle=\boldsymbol{Z}(s)+\boldsymbol{V}(s),$ (6b)

with the deterministic and stochastic transfer functions $G(s)$ and $H(s)$

$\displaystyle G(s)=\begin{bmatrix}g_{11}(s)&\cdots&g_{1n_{u}}(s)\\\
\vdots&\ddots&\vdots\\\
g_{n_{z}1}(s)&\cdots&g_{n_{z}n_{u}}(s)\end{bmatrix},\,g_{ij}(s)=\frac{b_{ij}(s)}{a_{ij}(s)}e^{-\tau_{ij}s},$
(7a) $\displaystyle H(s)=\begin{bmatrix}h_{11}(s)&\cdots&h_{1n_{w}}(s)\\\
\vdots&\ddots&\vdots\\\
h_{n_{z}1}(s)&\cdots&h_{n_{z}n_{w}}(s)\end{bmatrix},\,h_{ip}(s)=\frac{d_{ip}(s)}{c_{ip}(s)},$
(7b)

where $U$ and $\boldsymbol{W}$ are the deterministic and stochastic inputs to
the system. $\boldsymbol{Z}$ and $\boldsymbol{Y}$ are the output and
measurement and $\boldsymbol{V}$ indicates the measurement noise. We assume
that $g_{ij}(s)$ for $i=1,2,\ldots,n_{z}$, $j=1,2,\ldots,n_{u}$ with the input
time delay $\left\\{\tau_{ij}\right\\}\in\mathbb{R}_{0}^{+}$ is proper, and
$h_{ip}(s)$ for $p=1,2,\ldots,n_{w}$ is strictly proper and has no time delay.

Based on [8], using

$Z^{d}(s)=G(s)U(s),\qquad\boldsymbol{Z}^{s}(s)=H(s)\boldsymbol{W}(s),$ (8)

and $\boldsymbol{Z}(s)=Z^{d}(s)+\boldsymbol{Z}^{s}(s)$, the Noise-Separation
(NS) state space realization of (6) can be performed as

$\displaystyle Z^{d}(s)\sim\left\\{\begin{aligned}
\dot{x}^{d}_{ij}(t)&=A_{c,ij}^{d}x_{ij}^{d}(t)+B_{c,ij}^{d}u_{j}(t-\tau_{ij}),\\\
z^{d}_{ij}(t)&=C_{c,ij}^{d}x_{ij}^{d}(t)+D_{c,ij}^{d}u_{j}(t-\tau_{ij}),\end{aligned}\right.$
(9a) $\displaystyle\boldsymbol{Z}^{s}(s)\sim\left\\{\begin{aligned}
d\boldsymbol{x}^{s}(t)&=A_{c}^{s}\boldsymbol{x}^{s}(t)dt+B_{c}^{s}d\boldsymbol{\omega}(t),\\\
\boldsymbol{z}^{s}(t)&=C_{c}^{s}\boldsymbol{x}^{s}(t),\quad
d\boldsymbol{\omega}(t)\sim N_{iid}(0,Idt),\end{aligned}\right.$ (9b)

with the system variables

$\displaystyle
x^{d}=\begin{bmatrix}{x^{d}_{11}}^{\prime}&{x^{d}_{21}}^{\prime}&\ldots&{x^{d}_{n_{z}n_{u}}}^{\prime}\end{bmatrix}^{\prime},\quad
z^{d}_{i}(t)=\sum_{j=1}^{n_{u}}z^{d}_{ij}(t),$ (10a) $\displaystyle
z^{d}=\begin{bmatrix}z^{d}_{1}&z^{d}_{2}&\ldots&z^{d}_{n_{z}}\end{bmatrix}^{\prime},\qquad\boldsymbol{z}(t)=z^{d}(t)+\boldsymbol{z}^{s}(t),$
(10b)

where we convert the deterministic part $Z^{d}(s)$ into $[i\times j]$ SISO
time-delay state space models as they have different time delays and cannot be
expressed by one standard state space model. Note that
$\mathcal{G}_{c}=\left\\{A_{c,ij}^{d},B_{c,ij}^{d},C_{c,ij}^{d},D_{c,ij}^{d}\right\\}$
and $\mathcal{H}_{c}=\left\\{A_{c}^{s},B_{c}^{s},C_{c}^{s}\right\\}$ are
deterministic and stochastic system matrices.

The corresponding discrete-time state space models are

$\displaystyle Z^{d}(s)\sim\left\\{\begin{aligned}
&x^{d}_{k+1}=A^{d}x_{k}^{d}+B^{d}u_{k},\\\
&z^{d}_{k}=C^{d}x^{d}_{k}+D^{d}u_{k},\end{aligned}\right.$ (11a)
$\displaystyle\boldsymbol{Z}^{s}(s)\sim\left\\{\begin{aligned}
&\boldsymbol{x}^{s}_{k+1}=A^{s}\boldsymbol{x}_{k}^{s}+\boldsymbol{w}_{k},\\\
&\boldsymbol{z}^{s}_{k}=C^{s}\boldsymbol{x}^{s}_{k},\end{aligned}\right.$
(11b)

and

$\displaystyle\boldsymbol{y}_{k}=\boldsymbol{z}_{k}+\boldsymbol{v}_{k},\quad\begin{bmatrix}\boldsymbol{w}_{k}\\\
\boldsymbol{v}_{k}\end{bmatrix}\sim N_{iid}\left(\begin{bmatrix}0\\\
0\end{bmatrix},\begin{bmatrix}R_{ww}&0\\\ 0&R_{vv}\end{bmatrix}\right),$ (12)

where the deterministic system (11a) are obtained by stacking the all SISO
deterministic state space models [22]. Note that
$\mathcal{G}_{d}=\left\\{A^{d},B^{d},C^{d},D^{d}\right\\}$ and
$\mathcal{H}_{d}=\left\\{A^{s},C^{s}\right\\}$ are corresponding discrete
system matrices.

The covariance
$R_{ww}=\int_{0}^{T_{s}}e^{A_{c}^{s}t}B_{c}^{s}{B_{c}^{s}}^{\prime}{e^{A_{c}^{s}t}}^{\prime}dt$
can be solved using ordinary differential equation (ODE) methods or by
$\footnotesize{\begin{bmatrix}\Phi_{11}&\Phi_{12}\\\
0&\Phi_{21}\end{bmatrix}=exp(\begin{bmatrix}-A_{c}^{s}&B_{c}^{s}{B_{c}^{s}}^{\prime}\\\
0&{A_{c}^{s}}^{\prime}\end{bmatrix}})$ and
$R_{ww}=\Phi_{22}^{\prime}\Phi_{12}$.

The main advantage of using the NS state space realization is one may run the
filtering with (9b) and run the control algorithm with (9a). The Kalman filter
for (11b) can be performed as

$\displaystyle\hat{y}^{s}_{k|k-1}=C^{s}\hat{x}^{s}_{k|k-1},\qquad
e_{k}=y^{s}_{k}-\hat{y}^{s}_{k|k-1},$ (13a)
$\displaystyle\qquad\quad\hat{x}^{s}_{k|k}=\hat{x}^{s}_{k|k-1}+K_{fx}e_{k},$
(13b) with the measurement covariance and the Kalman gain
$R_{e}=C^{s}P({C^{s}})^{\prime}+R_{ww},\quad
K_{fx}=P({C^{s}})^{\prime}R_{e}^{-1},$ (13c)

where $P$ indicates the stationary state error covariance obtained by the
solution of discrete-time Algebraic Riccati Equation (DARE)

$P=A^{s}P({A^{s}})^{\prime}-A^{s}P({C^{s}})^{\prime}R_{e}^{-1}C^{s}P({A^{s}})^{\prime}+R_{vv}.$
(14)

We then could obtain the following prediction expressions

$\displaystyle\hat{x}^{d}_{k+j+1|k}=A^{d}\hat{x}^{d}_{k+j|k}+B^{d}u_{k+j},\>\>\>\>\hat{x}^{s}_{k+j+1|k}=A^{s}\hat{x}^{s}_{k+j|k},$
(15a)
$\displaystyle\hat{z}^{d}_{k+j|k}=C^{d}\hat{x}^{d}_{k+j|k}+D^{d}u_{k+j},\quad\>\>\>\hat{z}^{s}_{k+j|k}=C^{s}\hat{x}^{s}_{k+j|k},$
(15b)
$\displaystyle\hat{z}_{k+j|k}=\hat{z}^{d}_{k+j|k}+\hat{z}^{s}_{k+j|k},\qquad\qquad\text{for}\>\>j=1,\ldots,N.$
(15c)

######  Remark 1.

The NS state space expressions (9) corresponding to the linear stochastic
differential equation (SDE) model

$\displaystyle d\boldsymbol{x}(t)$
$\displaystyle=\left(A_{c}x(t)+B_{c}u(t-\tau)\right)dt+G_{c}d\boldsymbol{\omega}(t),$
(16a) $\displaystyle\boldsymbol{z}(t)$
$\displaystyle=C_{c}\boldsymbol{x}(t)+D_{c}u(t-\tau),$ (16b) where
$x=[x^{d};x^{s}]$ and $x\sim N(\hat{x}_{0},P_{0})$ and $\displaystyle A_{c}=$
$\displaystyle\begin{bmatrix}A^{d}_{c}&0\\\ 0&A^{s}_{c}\end{bmatrix},\quad
B_{c}=\begin{bmatrix}B^{d}_{c}\\\ 0\end{bmatrix},\quad
G_{c}=\begin{bmatrix}0\\\ B_{c}^{s}\end{bmatrix},$ (16c) $\displaystyle\quad
C=\begin{bmatrix}C^{d}&C^{s}\end{bmatrix},\qquad D_{c}=D^{d}_{c}.$ (16d)

Note that the above expressions are for the SISO case. The corresponding DT
stochastic state space model is

$\displaystyle\boldsymbol{x}_{k+1}$
$\displaystyle=A\boldsymbol{x}_{k}+Bu_{k}+G\boldsymbol{w}_{k},$ (17a)
$\displaystyle\boldsymbol{z}_{k}$ $\displaystyle=C\boldsymbol{x}_{k}+Du_{k}.$
(17b)

###  II-B Reference tracking and input regularization objectives

Define the output and input reference tracking error $\tilde{z}(t)$ as

$\tilde{z}(t)=\begin{bmatrix}z(t)-\bar{z}(t)\\\
u(t)-\bar{u}(t)\end{bmatrix}=\begin{bmatrix}z(t)-\bar{z}_{k}\\\
u(t)-\bar{u}_{k}\end{bmatrix},\quad\text{for}\>\>t_{k}\leq t<t_{k+1},$ (18)

where we consider piece-wise constant parameterization on references
$\bar{z}(t)=\bar{z}_{k}$ and $\bar{u}(t)=\bar{u}_{k}$ for $t_{k}\leq
t<t_{k+1}$.

We then can define the following CT LQ-OCP for minimizing the output and input
reference tracking error as

$\displaystyle\min_{x,u,z,\tilde{z}}\quad$
$\displaystyle\phi_{z}+\phi_{u}=\int_{t_{0}}^{t_{0}+T}l_{\tilde{z}}(\tilde{z}(t))dt$
(19a) $\displaystyle s.t.$
$\displaystyle\boldsymbol{x}(t_{0})=\hat{x}_{0},\leavevmode\nobreak\
\eqref{eq:ZOHInput},\leavevmode\nobreak\ \eqref{eq:ns-
DeterministicModel},\leavevmode\nobreak\ \eqref{eq:ns-
outputFunctions},\leavevmode\nobreak\ \eqref{eq:trackingError},$ (19b)

with the stage cost function $l_{\tilde{z}}(\tilde{z}(t))$

$l_{\tilde{z}}(\tilde{z}(t))=\frac{1}{2}\left\|W_{\tilde{z}}\tilde{z}(t)\right\|^{2}_{2}=\frac{1}{2}\tilde{z}(t)^{\prime}Q_{c\tilde{z}}\tilde{z}(t),$
(20)

where
$Q_{c\tilde{z}}=\text{diag}\left(Q_{cz},Q_{cu}\right)=W_{\tilde{z}}^{\prime}W_{\tilde{z}}$.
$Q_{cz}$ and $Q_{cu}$ are weight matrices for the reference tracking and input
regularization objectives.

The above problem is in CT and subject to a linear deterministic time-delay
system. The corresponding DT equivalent of (19) is

$\displaystyle\min_{x,u}\quad$
$\displaystyle\phi_{z}+\phi_{u}=\sum_{k\in\mathcal{N}}l_{\tilde{z},k}(x_{k},u_{k})$
(21a) $\displaystyle s.t.$ $\displaystyle
x_{0}=\hat{x}_{0},\leavevmode\nobreak\ \eqref{eq:ns-discreteStateSpace-
determinsitic},\leavevmode\nobreak\ \eqref{eq:ns-predictionProcess},$ (21b)

with the stage cost function $l_{\tilde{z},k}(x_{k},u_{k})$

$l_{\bar{z},k}(x_{k},u_{k})=\frac{1}{2}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}^{\prime}Q\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}+q_{k}^{\prime}\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix}+\rho_{k},$ (22)

where

$q_{k}=M\left[\bar{z}_{k};\bar{u}_{k}\right],\quad\rho_{k}=\frac{1}{2}\begin{bmatrix}\bar{z}_{k}\\\
\bar{u}_{k}\end{bmatrix}^{\prime}Q_{c\tilde{z}}\begin{bmatrix}\bar{z}_{k}\\\
\bar{u}_{k}\end{bmatrix}T_{s},\quad k\in\mathcal{N}.$ (23)

Note that when using NS state space realization, the reference trajectory
becomes $\bar{z}_{k+j}=\bar{z}(t)-\hat{z}^{s}_{k+j|k}$ for
$t\in[t_{k+j},t_{k+j+1})$.

######  Remark 2.

The studies [23, 22] described the discrete system matrices $(A,B,Q,M)$ as
differential equations and proposed a novel step-doubling method for solving
them. Their numerical experiments showed that the proposed numerical method is
significantly faster than the standard ODE method while keeping the same
accuracy.

Note that $u=[u_{0};u_{1};\cdots;u_{N-1}]$ and $u_{k}=I_{k}u$ for
$I_{k}=[0\cdots I\cdots 0]$ and $k\in\mathcal{N}$. The system state $x_{k}$
can be expressed as

$x_{k}=b_{k}+\Gamma_{k}u,\quad
b_{k}=A^{k}x_{0},\quad\Gamma_{k}=\sum_{i=0}^{k}A^{k-1-i}BI_{i}.$ (24)

We then rewrite the LQ-OCP (21) in the form of QP

$\displaystyle\phi_{z}+\phi_{u}=\frac{1}{2}u^{\prime}H_{\tilde{z}}u+g_{\tilde{z}}^{\prime}u,$
(25)

where

$H_{\tilde{z}}=\sum_{k=0}^{N-1}\begin{bmatrix}\Gamma_{k}\\\
I_{k}\end{bmatrix}^{\prime}Q\begin{bmatrix}\Gamma_{k}\\\
I_{k}\end{bmatrix},\>\>g_{\tilde{z}}=\sum_{k=0}^{N-1}\begin{bmatrix}\Gamma_{k}\\\
I_{k}\end{bmatrix}^{\prime}\left(Q\begin{bmatrix}b_{k}\\\
0\end{bmatrix}+q_{k}\right).$ (26)

###  II-C Input ROM and economic objectives

We then consider the LQ-OCP for penalizing the input rate-of-movement (ROM)
and input cost (economic)

$\begin{split}\phi_{\Delta
u}+\phi_{eco}&=\int_{t_{0}}^{t_{0}+T}l_{\tilde{u}}(\dot{u}(t),u(t))dt,\end{split}$
(27)

and the stage cost function is

$l_{\tilde{u}}(\dot{u}(t),u(t))=\frac{1}{2}\left\|W_{c\Delta
u}\dot{u}(t)\right\|^{2}_{2}+q_{ceco}^{\prime}u(t),$ (28)

where $Q_{c\Delta u}=W_{c\Delta u}^{\prime}W_{c\Delta u}$ and $q_{ceco}$ are
the weight matrices of the input ROM and the economic objectives.

The discretization of the input cost objective is intuitive if a zero-order
hold (ZOH) parameterization is applied to the inputs. However, this is not the
case for the input ROM objective function, as the input derivative is not even
defined. The previous work [7] described the discretization schemes on the
input ROM penalty function using piece-wise affine functions and the ZOH
discretization. The author approved that the solution of the discrete problem
is convergent in $L^{2}$ towards the original problem’s minimizers. When
taking the number of the discretization $N\rightarrow\infty$, the minimizers
of the discrete problem converge in $H^{1}$ towards the solution of the
original problem.

Lets consider the ZOH discretization on the input $u(t)=u_{k}$ for
$t\in[t_{k},t_{k+1})$ and the discrete approximation of the input ROM penalty
is

$\phi_{\Delta
u}=\frac{1}{T_{s}}\sum_{k\in\mathcal{N}}\left\|u_{k}-u_{k-1}\right\|^{2}_{Q_{c\Delta
u}}.$ (29)

We then can perform the discretization of (27) as

$\phi_{\Delta
u}+\phi_{eco}=\sum_{k\in\mathcal{N}}l_{\tilde{u},k}(u_{k},u_{k-1}),$ (30)

with the stage cost function

$\begin{split}l_{\tilde{u},k}(u_{k},u_{k-1})=\frac{1}{2}\begin{bmatrix}u_{k}\\\
u_{k-1}\end{bmatrix}^{\prime}\bar{Q}_{\Delta u}\begin{bmatrix}u_{k}\\\
u_{k-1}\end{bmatrix}+q_{eco}^{\prime}u_{k},\end{split}$ (31)

where $u_{-1}$ indicates the input at $t_{0-1}$. The discrete-time weights
$\bar{Q}_{\Delta u}$ and $q_{eco}$ are

$\bar{Q}_{\Delta u}=\begin{bmatrix}Q_{\Delta u}&-Q_{\Delta u}\\\ -Q_{\Delta
u}&Q_{\Delta u}\end{bmatrix},\>Q_{\Delta u}=\frac{Q_{c\Delta
u}}{T_{s}},\>q_{eco}=q_{ceco}T_{s}.$ (32)

Consequently, along with possible input box and input ROM constraints, we can
express the input ROM and economic objective functions in the form of QP as

$\displaystyle\min_{u}\quad$ $\displaystyle\phi_{\Delta
u}+\phi_{u}=\frac{1}{2}u^{\prime}H_{\tilde{u}}u+g_{\tilde{u}}^{\prime}u$ (33a)
$\displaystyle s.t.$ $\displaystyle u_{\min,k}\leq u_{k}\leq u_{\max,k},\quad$
$\displaystyle k\in\mathcal{N},$ (33b) $\displaystyle\Delta
u_{\min,k}\leq\Delta u_{k}\leq\Delta u_{\max,k},\quad$ $\displaystyle
k\in\mathcal{N},$ (33c)

where the quadratic and linear term coefficients

$\displaystyle H_{\tilde{u}}=\sum_{k=1}^{N-1}\begin{bmatrix}I_{k}\\\
I_{k-1}\end{bmatrix}^{\prime}\bar{Q}_{\Delta u}\begin{bmatrix}I_{k}\\\
I_{k-1}\end{bmatrix}+I_{0}^{\prime}Q_{\Delta u}I_{0},$ (34a) $\displaystyle
g_{\tilde{u}}=\sum_{k=0}^{N-1}-I_{0}^{\prime}Q_{\Delta
u}u_{-1}+I_{k}^{\prime}q_{eco}.$ (34b)

###  II-D Soft output constraint penalty

We then introduce the soft output constraints

$\displaystyle z_{k+j|k}\geq z_{\min,k+j|k}-\xi_{k+j},$ $\displaystyle
k=1,2,\ldots,N,$ (35a) $\displaystyle z_{k+j|k}\leq
z_{\max,k+j|k}+\eta_{k+j},$ $\displaystyle k=1,2,\ldots,N,$ (35b)
$\displaystyle\xi_{k+j}\geq 0,$ $\displaystyle k=1,2,\ldots,N,$ (35c)
$\displaystyle\eta_{k+j}\geq 0,$ $\displaystyle k=1,2,\ldots,N,$ (35d) where
$\xi$ and $\eta$ are slackness variables. The output $z_{k}$ is subject to the
deterministic system (11). $z_{\min,k+j|k}=z_{\min,k+j}-\hat{z}^{s}_{k+j|k}$
and $z_{\max,k+j|k}=z_{\max,k+j}-\hat{z}^{s}_{k+j|k}$ are modified soft
constraints.

The corresponding penalty function

$\displaystyle\phi_{\xi}+\phi_{\eta}=\int_{t_{0}}^{t_{0}+T}l_{\xi}(\xi(t))+l_{\eta}(\eta(t))dt,$
(36)

with the stage cost functions

$\displaystyle
l_{\xi}(\xi(t))=\frac{1}{2}\left\|W_{c\xi}\xi(t)\right\|^{2}_{2}+q_{c\xi}^{\prime}\xi(t),$
(37) $\displaystyle
l_{\eta}(\eta(t))=\frac{1}{2}\left\|W_{c\eta}\eta(t)\right\|^{2}_{2}+q_{c\eta}^{\prime}\eta(t),$
(38)

where $Q_{c\xi}=W_{c\xi}^{\prime}W_{c\xi}$,
$Q_{c\eta}=W_{c\eta}^{\prime}W_{c\eta}$, $q_{c\xi}$ and $q_{c\eta}$ are weight
matrices of the quadratic and linear penalty functions, respectively. We
assume piece-wise constant $\xi(t)=\xi_{k}$ and $\eta(t)=\eta_{k}$ for
$t\in[t_{k},\,t_{k}+T_{s})$. The corresponding discrete equivalent is

$\displaystyle\phi_{\xi}+\phi_{\eta}=\sum_{k=1}^{N}\frac{1}{2}(\left\|\xi_{k}\right\|^{2}_{Q_{\xi}}+\left\|\eta_{k}\right\|^{2}_{Q_{\eta}})+q_{\xi}^{\prime}\xi_{k}+q_{\eta}^{\prime}\eta_{k},$
(39)

where $Q_{\xi}=T_{s}Q_{c\xi}$, $q_{\xi}=T_{s}q_{c\xi}$,
$Q_{\eta}=T_{s}Q_{c\eta}$ and $q_{\eta}=T_{s}q_{c\eta}$.

The DT penalty functions (39) can be rewritten in the form of QP as

$\displaystyle\phi_{\xi}+\phi_{\eta}=\frac{1}{2}\begin{bmatrix}\xi\\\
\eta\end{bmatrix}^{\prime}H_{\tilde{s}}\begin{bmatrix}\xi\\\
\eta\end{bmatrix}+g_{\tilde{s}}^{\prime}\begin{bmatrix}\xi\\\
\eta\end{bmatrix},$ (40a)

where $\xi$=$[\xi_{1};\xi_{2};\cdots;\xi_{N}]$ and
$\eta$=$[\eta_{1};\eta_{2};\cdots;\eta_{N}]$ and

$\displaystyle
H_{\xi}=\text{diag}({Q}_{\xi};{Q}_{\xi};\cdots;Q_{\xi}),\>g_{\xi}=\begin{bmatrix}q_{\xi}&q_{\xi}&\cdots&q_{\xi}\end{bmatrix}^{\prime},$
(41a) $\displaystyle
H_{\eta}=\text{diag}({Q}_{\eta};{Q}_{\eta};\cdots;Q_{\eta}),g_{\eta}=\begin{bmatrix}q_{\eta}&q_{\eta}&\cdots&q_{\eta}\end{bmatrix}^{\prime},$
(41b) $\displaystyle
H_{\tilde{s}}=\text{diag}({H}_{\xi};{H}_{\eta}),\qquad\quad\>\>\>g_{\tilde{s}}=\begin{bmatrix}g_{\xi}&g_{\eta}\end{bmatrix}^{\prime}.$
(41c)

###  II-E Design and implementation of CT-LMPC

Algorithm 1 Design of CT-LMPC

Require: $(G(s),H(s),Q_{c\tilde{z}},Q_{c\Delta
u},Q_{c\xi},Q_{c\eta},q_{ceco},q_{c\xi},q_{c\eta},$

$R_{vv},N,T_{s})$

Compute NS state space matrices $\mathcal{G}_{c}$, $\mathcal{H}_{c}$,
$\mathcal{G}_{d}$ and $\mathcal{H}_{d}$ using (9) and (11)

Compute Kalman gain $K_{fx}$ using (14)

Compute DT weight matrices $Q$ and $M$ using (22)

Compute $b=\left\\{b_{k}\right\\}$ and $\Gamma=\left\\{\Gamma_{k}\right\\}$
for $k\in\mathcal{N}$ using (24)

Compute $H_{\tilde{z}}$ of $\phi_{z}$ and $\phi_{u}$ using (26)

Compute $H_{\tilde{u}}$ of $\phi_{\Delta u}$ and $\phi_{eco}$ using (34a)

Compute $H_{\tilde{s}}$ and $g_{\tilde{s}}$ of $\phi_{\xi}$ and $\phi_{\eta}$
using (40)

Set $H=$ diag($[H_{\tilde{z}}+H_{\tilde{u}};H_{\tilde{s}}]$)

Return:
$(\mathcal{G}_{d},\mathcal{H}_{d},K_{fx},Q,M,b,\Gamma,H,g_{\tilde{s}})$

Algorithm 2 Implementation of CT-LMPC

Require: $(\bar{z},\bar{u},u_{\min},u_{\max},\Delta u_{\min},\Delta
u_{max},z_{\min},z_{\max}$,

$u_{k-1},y_{k},\hat{x}^{s}_{k-1|k-1},\mathcal{G}_{d},\mathcal{H}_{d},K_{fx},Q,M,b,\Gamma,H,g_{\tilde{s}})$

Run Kalman filter for $\hat{x}^{s}_{k|k}$ and
$\left\\{\hat{z}^{s}_{k+j|k}\right\\}$ using (15b)

for j=1,2,$\ldots$,N do

Update $\bar{z}_{k+j}$, $z_{\min,k+j}$ and $z_{\max,k+j}$ with
$\hat{z}^{s}_{k+j|k}$

Set $q_{k+j}=M[\bar{z}_{k+j};\bar{u}_{k+j}]$

end for

Compute $g_{\tilde{z}}$ using (26)

Compute $g_{\tilde{u}}$ using (34b)

Set $g=[g_{\tilde{z}}+g_{\tilde{u}};g_{\tilde{s}}]$

Solve the QP (42) to get the optimal solution $u^{*}$

Let $u_{k}\leftarrow u^{*}(1:n_{u})$

Return: $(u_{k},\hat{x}^{s}_{k|k})$

Combining the objective functions introduced previously, we have

$\displaystyle\min_{\left\\{u,\xi,\eta\right\\}}\quad$
$\displaystyle\phi=\phi_{z}+\phi_{u}+\phi_{\Delta
u}+\phi_{eco}+\phi_{\xi}+\phi_{\eta}$ (42a) $\displaystyle s.t.$
$\displaystyle x_{0}=\hat{x}_{0},\leavevmode\nobreak\ \eqref{eq:ns-
discreteStateSpace-determinsitic},\leavevmode\nobreak\
\eqref{eq:inputConstraints},\leavevmode\nobreak\
\eqref{eq:inputROMConstraints},\leavevmode\nobreak\
\eqref{eq:softConstraints},$ (42b)

where the objective functions $\phi_{z}$, $\phi_{u}$, $\phi_{\Delta u}$,
$\phi_{eco}$, $\phi_{\xi}$ and $\phi_{\eta}$ are the corresponding discrete
equivalent of their original CT LQ-OCPs.

The cost function $\phi$ can be expressed in the form of QP as

$\phi=\frac{1}{2}\begin{bmatrix}u\\\ \xi\\\
\eta\end{bmatrix}^{\prime}H\begin{bmatrix}u\\\ \xi\\\
\eta\end{bmatrix}+g^{\prime}\begin{bmatrix}u\\\ \xi\\\ \eta\end{bmatrix},$
(43)

where

$H=\begin{bmatrix}H_{\tilde{z}}+H_{\tilde{u}}&0\\\
0&H_{\tilde{s}}\end{bmatrix},\qquad
g=\begin{bmatrix}g_{\tilde{z}}+g_{\tilde{u}}\\\ g_{\tilde{s}}\end{bmatrix}.$
(44)

Consequently, we obtain the objective function (42) that is the discrete-time
equivalent of the CT LQ-OCPs introduced in previous subsections. Algorithm 1
describes the design of CT-LMPC and the implementation of the CT-LMPC is
illustrated by Algorithm 2.

(a) Deterministic sim. $T_{s}$ = 5 [s]

(b) Deterministic sim. $T_{s}$ = 15 [s]

(c) Deterministic sim. $T_{s}$ = 25 [s]

(d) Stochastic sim. $T_{s}$ = 5 [s]

(e) Stochastic sim. $T_{s}$ = 15 [s]

(f) Stochastic sim. $T_{s}$ = 25 [s]

Figure 1: SISO example closed-loop simulations with different controller
sampling times $T^{c}_{s}=5,15,25$ [s]. The blue curves indicate the results
obtained from DT-LMPC and the black curves are the results of CT-LMPC.

##  III Numerical Experiments

This section describes closed-loop simulations with the proposed CT-LMPC and
standard DT-LMPC. Consider the simulation model with the following expressions

$\displaystyle\boldsymbol{Z}(s)=G(s)U(s)+G_{d}(s)(D(s)+\boldsymbol{W}(s)),$
(45a) $\displaystyle\boldsymbol{Y}(s)=\boldsymbol{Z}(s)+\boldsymbol{V}(s),$
(45b)

and it can be converted to a state space model with the sampling time $T_{s}$
as

$\displaystyle\boldsymbol{x}_{k+1}$
$\displaystyle=A\boldsymbol{x}_{k}+Bu_{k}+Ed_{k}+G\boldsymbol{w}_{k},\>\>\boldsymbol{w}_{k}\sim
N_{iid}(0,R_{ww}),$ (46a) $\displaystyle\boldsymbol{y}_{k}$
$\displaystyle=C\boldsymbol{x}_{k}+Du_{k}+\boldsymbol{v}_{k},\qquad\qquad\>\>\boldsymbol{v}_{k}\sim
N_{iid}(0,R_{vv}),$ (46b)

where $x_{k}\sim N(\hat{x}_{k},P)$ is the states, $u_{k}$ is the inputs and
$d_{k}$ is the disturbance.

In the closed-loop simulations, we use (46) as the simulator. We develop the
DT-LMPC based on the previous work by [8]. We estimate the control model (6)
for CT- and DT-LMPC by step-response experiments. It may be done using system
parameters identification schemes introduced by [10, 11, 15].

###  III-A Closed-loop simulation - a SISO example

In the following, we perform a series of closed-loop simulations with CT- and
DT-LMPC implementations for a SISO system. The SISO example considers the
simulation model (45) with transfer functions

$\displaystyle g(s)=\frac{10.12(-3.41s+1)}{(15.9s+1)(24.2s+1)},$ (47a)
$\displaystyle g_{d}(s)=\frac{-0.5}{(5.8s+1)(4.7s+1)},$ (47b)

and we convert the transfer functions into the DT state space (46) with
sampling time $T_{s}=1$ [s]. The unknown disturbance $d_{k}$ = 2.0 for
$3.75\leq t\leq 11.25$ [min]. The system covariances are
$R_{ww}=R_{vv}=0.05^{2}$.

The estimated control model (6) for the two MPCs is

$\hat{g}(s)=\frac{10.12(-3.58s+1)}{(18.9s+1)(22.2s+1)},\>\>\hat{h}(s)=\frac{1}{s}\frac{0.6}{(s+1)},$
(48)

where we select $\hat{h}(s)$ as an integrator to ensure the controller can
handle step-type error [14, 18].

In the SISO example, we consider the objective function (42) with only the
reference tracking $\phi_{z}$ and input ROM $\phi_{\Delta u}$ cost functions.
For DT-LMPC, $\phi_{z}$ and $\phi_{\Delta u}$ are designed as

$\phi=\sum_{k=0}^{N-1}\left\|z_{k+1}-\bar{z}_{k+1}\right\|^{2}_{Q_{cz}}+\left\|u_{k}-u_{k-1}\right\|^{2}_{Q_{c\Delta
u}},$ (49)

The prediction and control horizon is $N=20$. The output target is $\bar{z}=2$
for $t\leq 7.5$ [min]. The weight matrices of the two objectives are
$Q_{cz}=20$ and $Q_{c\Delta u}=1$, and the input constraints are $u_{\min}$ =
-1 and $u_{\max}$ = 1.

(a) The outputs of deterministic and stochastic sim.

(b) The inputs of deterministic and stochastic sim.

Figure 2: The closed-loop simulations of a simulated cement mill system with
DT-LMPC (blue curves) and CT-LMPC (black curves) implementations. The
simulator and controller sampling times are $T_{s}$ = 1 [min] and $T_{s}^{c}$
= 2 [min]

To demonstrate the contrast between CT- and DT-LMPC, we discretize both using
different controller sampling times $T^{c}_{s}$ = 5, 15, 25 [s]. Fig. 1 shows
closed-loop simulations of the SISO system with CT- and DT-LMPC
implementations. In Fig. 1(a) and 1(d), both CT-LMPC (black curves) and DT-
LMPC (blue curves) effectively control the SISO system, allowing the outputs
to track the given target (red dashed line). Two overshoots occur at $t=3.75$
and $11.25$ [min] due to an unknown disturbance, and the two MPCs can handle
it after a few iterations. In this case, the simulation results show little
difference between the two MPCs. Upon increasing $T^{c}_{s}$ to 15 [s],
additional overshoots emerge around 1 and 7.5 [min] in Fig. 1(b) and 1(e).
Notably, the second overshoot at $t=7.5$ [min] for DT-LMPC is larger than that
for CT-LMPC. Moreover, the overshoots resulting from the unknown disturbance
are more significant compared to the previous scenario, yet both MPCs can
manage it. Fig. 1(c) and 1(f) show the closed-loop simulations when
$T_{c}^{s}$ = 25 [s]. Here, we can see a significant decrease in the closed-
loop performance of both MPCs compared to the previous two cases. Both CT- and
DT-LMPC exhibit oscillations and struggle to follow the target. However, the
amplitude of the output from CT-LMPC is considerably smaller than that from
DT-LMPC. CT-LMPC maintains control of the output to the desired target after
1-2 iterations, while DT-LMPC fails to do so.

The SISO example simulation results show that a large sampling time can
degrade the closed-loop performance of the predictive controller and even lead
to system instability. For the same system parameters, the proposed CT-LMPC
has a better closed-loop performance than the traditional DT-LMPC, and this
gap will increase with the increase of the sampling time.

###  III-B Closed loop simulation - a MIMO example

The MIMO example concerns the cement mill system introduced in [19, 15], where
the cement mill system is described by transfer functions

$\displaystyle
G(s)=\begin{bmatrix}\frac{0.62e^{-5s}}{(45s+1)(8s+1)}&\frac{0.29(8s+1)e^{-1.5s}}{(2s+1)(38s+1)}\\\
\frac{-15e^{-5s}}{60s+1}&\frac{5e^{-0.1s}}{(14s+1)(s+1)}\end{bmatrix},$ (50a)
$\displaystyle G_{d}(s)=\begin{bmatrix}\frac{-1.0e^{-3s}}{(32s+1)(21s+1)}\\\
\frac{60}{(30s+1)(20s+1)}\end{bmatrix}.$ (50b)

The system inputs $u$ = [feed flow rate (TPH); separator speed (%)] and the
system outputs $z$ = [elevator load (kW); fineness (cm2/g)]. The system
disturbance $D$ is the clinker hardness (HGI). The transfer functions of the
estimated control model (6) for the two MPCs are

$\displaystyle\hat{G}(s)=\begin{bmatrix}\frac{0.8e^{-4s}}{(30s+1)(15s+1)}&\frac{-0.45e^{-2s}}{30.0s+1}\\\
\frac{-17.7e^{-5s}}{(65s+1)(15s+1)}&\frac{-19.4}{15s+1}\end{bmatrix},$ (51a)
$\displaystyle\hat{H}(s)=\text{diag}\left(\left[\frac{1}{s}\frac{0.5}{s+1};\frac{1}{s}\frac{1}{s+1}\right]\right),$
(51b)

where we manually introduce the plant-model mismatch on the deterministic
model $\hat{G}(s)$.

The cement mill system is converted into a DT state space model (46) with the
sampling time $T_{s}$ = 1 [min]. The simulation time $T_{sim}=1.5$ [h]. The
related covariances are selected as $R_{ww}$ = 1.0 and $R_{vv}$ = diag([0.1;
50]). The unknown disturbance $d_{k}$ = 8 for $0$[min] and $d_{k}$ = 0 for the
rest of the time.

For the CT- and DT-LMPC, the controller sampling time $T_{s}^{c}=2$ [min], the
prediction and control horizon is $N=60$. The weight matrices for the
reference tracking, input ROM, economics and soft output constraints penalty
are $Q_{cz}=$diag($200,10$), $Q_{c\Delta u}=$ diag($2,1$), $q_{ceco}=[2;1]$,
$Q_{c\xi}=Q_{c\eta}=$ diag($2000,100$) and $q_{c\xi}=q_{c\eta}=[20;1]$,
respectively. The input hard constraints are $u_{\min}=[-10;-20]$,
$u_{\max}=[10,20]$, $\Delta u_{\min}=[-5,-10]$ and $\Delta u_{\max}=[5,10]$.
The soft output constraints $z_{\min}=[-2,-20]$ and $z_{\max}=[2,20]$ and are
deviation variables.

Fig. 2 illustrates the closed-loop simulations of a simulated cement mill
system with CT- and DT-LMPC implementations. In Fig. 2(a), both CT-LMPC (black
curves) and DT-LMPC (blue curves) can regulate the outputs towards the desired
targets (blue dashed lines) and within the given bounds for most of the time.
An overshoot, attributed to the plant-model mismatch in $y_{2}$ (Fineness), is
observed. The two predictive controllers can handle the problem and bring the
output back to the target after a few iterations. At $t=0.5$ and $1.5$ hours,
two overshoots on the output occur due to the unknown disturbance. Both CT-
and DT-LMPC can reject the unknown disturbance. Additionally, an overshoot
arises at around t = 0.75 [h] on the Fineness due to a step change in the
output targets. DT-LMPC exhibits a slightly larger overshoot than CT-LMPC.

The MIMO example closed-loop simulations show that the proposed CT-LMPC and
DT-LMPC exhibit highly similar closed-loop performance. They even converge to
identical optimal solutions when using appropriate sampling time and tuning
parameters. This outcome aligns with the earlier SISO example presented in
Fig. 1(a) and 1(d), reinforcing the consistency of our findings.

##  IV Conclusions

This paper presents the design, discretization and implementation of CT-LMPC.
We introduce different objective functions of CT-LMPC and show their
discretization and implementation. From the numerical experiments, we notice
that the proposed CT-LMPC has the following features:

  * 1.

When using an appropriate sampling time and well-tuned weight matrices, the
proposed CT-LMPC has a highly similar control performance as the traditional
DT-LMPC, and they will converge to an identical (very similar) optimal
solution

  * 2.

The proposed CT-LMPC has a better closed-loop performance than the traditional
DT-LMPC when the controller sampling time is large. This gap will increase
with the increase in controller sampling time

## References

  * [1] W. Alt, C. Schneider, and M. Seydenschwanz, “Regularization and implicit euler discretization of linear-quadratic optimal control problems with bang-bang solutions,” _Applied Mathematics and Computation_ , vol. 287-288, pp. 104–124, 2016.
  * [2] K. J. Åström, _Introduction to stochastic control theory_. Academic Press New York, 1970.
  * [3] K. J. Åström and B. Wittenmark, _Computer-Controlled Systems: Theory and Design, Third Edition_ , ser. Dover Books on Electrical Engineering. Dover Publications, 2011.
  * [4] M. Cannon and B. Kouvaritakis, “Infinite horizon predictive control of constrained continuous-time linear systems,” _Automatica_ , vol. 36, no. 7, pp. 943–955, 2000.
  * [5] G. F. Franklin, J. D. Powell, and M. L. Workman, _Digital Control of Dynamic Systems_ , 2nd ed., ser. Electrical and Computer Engineering; Control Engineering. Reading Massachusetts: Addison-Wesley, 1990.
  * [6] G. Frison and J. B. Jørgensen, “Efficient implementation of the riccati recursion for solving linear-quadratic control problems,” in _2013 IEEE International Conference on Control Applications (CCA)_ , 2013, pp. 1117–1122.
  * [7] M. Hagdrup, “Model predictive control for systems described by stochastic diﬀerential equations,” Ph.D. dissertation, Technical University of Denmark, Kgs. Lyngby, Denmark, 2019.
  * [8] M. Hagdrup, D. Boiroux, Z. Mahmoudi, H. Madsen, N. Kjølstad Poulsen, B. Poulsen, and J. B. Jørgensen, “On the significance of the noise model for the performance of a linear MPC in closed-loop operation,” _IFAC-PapersOnLine_ , vol. 49, no. 7, pp. 171–176, 2016, 11th IFAC Symposium on Dynamics and Control of Process SystemsIncluding Biosystems DYCOPS-CAB 2016.
  * [9] W. W. Hager, “Runge-kutta methods in optimal control and the transformed adjoint system,” _Numerische Mathematik_ , vol. 87, pp. 247–282, 2000.
  * [10] J. B. Jørgensen and S. B. Jørgensen, “Comparison of Prediction-Error-Modelling Criteria,” in _2007 American Control Conference_ , 2007, pp. 140–146.
  * [11] ——, “MPC-Relevant Prediction-Error Identification,” in _2007 American Control Conference_ , 2007, pp. 128–133.
  * [12] H. Li and Y. Shi, “Event-triggered robust model predictive control of continuous-time nonlinear systems,” _Automatica_ , vol. 50, no. 5, pp. 1507–1513, 2014.
  * [13] M. M. Morato, J. E. Normey-Rico, and O. Sename, “Model predictive control design for linear parameter varying systems: A survey,” _Annual Reviews in Control_ , vol. 49, pp. 64–80, 2020.
  * [14] K. R. Muske and T. A. Badgwell, “Disturbance modeling for offset-free linear model predictive control,” _Journal of Process Control_ , vol. 12, no. 5, pp. 617–632, 2002.
  * [15] D. H. Olesen, J. K. Huusom, and J. B. Jørgensen, “A tuning procedure for ARX-based MPC of multivariate processes,” in _2013 American Control Conference_ , 2013, pp. 1721–1726.
  * [16] G. Pannocchia, J. B. Rawlings, D. Q. Mayne, and G. M. Mancuso, “Whither discrete time model predictive control?” _IEEE Transactions on Automatic Control_ , vol. 60, no. 1, pp. 246–252, 2015.
  * [17] G. Pannocchia, J. B. Rawlings, D. Q. Mayne, and W. Marquardt, “On computing solutions to the continuous time constrained linear quadratic regulator,” _IEEE Transactions on Automatic Control_ , vol. 55, no. 9, pp. 2192–2198, 2010.
  * [18] G. Pannocchia and J. Rawlings, “Disturbance models for offset-free model-predictive control,” _AIChE Journal_ , vol. 49, pp. 426–437, 02 2003.
  * [19] G. Prasath, B. Recke, M. Chidambaram, and J. B. Jørgensen, “Application of soft constrained MPC to a cement mill circuit,” _IFAC Proceedings Volumes_ , vol. 43, no. 5, pp. 302–307, 2010, 9th IFAC Symposium on Dynamics and Control of Process Systems.
  * [20] V. Sakizlis, J. Perkins, and E. Pistikopoulos, “Explicit solutions to optimal control problems for constrained continuous-time linear systems,” _IEE Proceedings-Control Theory and Applications_ , vol. 152, no. 4, pp. 443–452, 2005.
  * [21] L. Wang, “Continuous time model predictive control design using orthonormal functions,” _International Journal of Control_ , vol. 74, no. 16, pp. 1588–1600, 2001.
  * [22] Z. Zhang, S. Hørsholt, and J. B. Jørgensen, “Numerical discretization methods for linear quadratic control problems with time delays,” in _12th IFAC Symposium on Advanced Control of Chemical Proccesses (ADCHEM 2024) (accepted)_ , 2023.
  * [23] Z. Zhang, J. L. Svensen, M. W. Kaysfeld, A. H. D. Andersen, S. Hørsholt, and J. B. Jørgensen, “Numerical discretization methods for the extended linear quadratic control problem,” in _European Control Conference (ECC) 2024 (accepted)_ , 2023.