# Numerical Discretization Methods
for the Extended Linear Quadratic Control Problem

Zhanhao Zhang, Jan Lorenz Svensen, Morten Wahlgreen Kaysfeld,
Anders Hilmar Damm Christensen, Steen Hørsholt, John Bagterp Jørgensen

###### Abstract

In this study, we introduce numerical methods for discretizing continuous-time
linear-quadratic optimal control problems (LQ-OCPs). The discretization of
continuous-time LQ-OCPs is formulated into differential equation systems, and
we can obtain the discrete equivalent by solving these systems. We present the
ordinary differential equation (ODE), matrix exponential, and a novel step-
doubling method for the discretization of LQ-OCPs. Utilizing Euler-Maruyama
discretization with a fine step, we reformulate the costs of continuous-time
stochastic LQ-OCPs into a quadratic form, and show that the stochastic cost
follows the $\chi^{2}$ distribution. In the numerical experiment, we test and
compare the proposed numerical methods. The results ensure that the discrete-
time LQ-OCP derived using the proposed numerical methods is equivalent to the
original problem.

##  I Introduction

In optimal control theory, linear-quadratic (LQ) optimal control problems are
considered the fundamental but core problem. They involve a quadratic cost
function that needs to be optimized, and the controlled system is linear. Due
to the simplicity and analytical solvability of LQ-OCPs, it has widespread
practical applications in numerous fields, including engineering, aerospace,
biology and economics [5, 7, 24, 25]. Implementing continuous-time LQ-OCPs in
real-world scenarios might be infeasible, primarily because many practical
applications operate on digital platforms in digital form. Consequently,
discretization techniques become imperative in practice.

There is rich research on discretization and numerical solution methods for
optimal control problems [8, 13, 20, 10, 15, 26, 35]. However, this is not the
case for how to obtain the discrete-time equivalent from a continuous-time LQ-
OCP. Based on existing literature, there are two popular ways for LQ
discretization: 1) design a continuous-time LQ-OCP and subsequently derive its
discrete-time equivalent, 2) initially discretizing the continuous-time
system, followed by designing a discrete-time LQ-OCP with the discrete-time
system. Compared with the latter, the first method provides a better
approximation, which eventually leads to the solution of the original problem.

The continuous-time LQ-OCP can be converted into discrete-time by assuming
zero-order hold (ZOH) on the state vector. However, this approximation is a
crude approximation that is often inaccurate for large sample times [16]. In
[3, 4], the cost equivalent is obtained by extending the continuous-time LQ-
OCP’s cost function. It leads to the analytic expressions of the desired
equivalent, discrete weighting matrices. Further, a matrix exponential method
is introduced for calculating the cost equivalent [11, 1, 2, 29, 34].
Modelling sampled-data systems with traditional approaches has fundamental
difficulties, which can be resolved using incremental models. Incremental
models provide a seamless connection between continuous- and discrete-time
systems, and they can be implemented for optimal filtering and control [36,
14, 27, 31]. In addition, LQ-OCPs are associated with advanced control
algorithms, such as model predictive control (MPC). Prediction-Error-Methods
for identification of these models exist [21, 23, 22, 19]. Also efficient
numerical methods for the solution of extended LQ-OCPs exists, e.g. structure
exploiting algorithms that use a Riccati recursion [20, 12, 19]. LQ-OCPs
become hard to solve and analyze when the controlled systems are stochastic.
The analytic expression of the stochastic cost and its expectation can be
calculated using It$\hat{o}$ calculus [3, 4]. Many stochastic LQ-OCPs aim to
optimize the mean value of their costs, which may not be true for some
scenarios, e.g., conditional Value-at-Risk optimization problems (CVaR
optimization) [32, 30, 33, 6, 18]. Therefore, it is critical to investigate
the cost function distribution of stochastic LQ-OCPs.

The key problem that we address in this paper:

  * 1.

Formulation of differential equation systems for LQ discretization

  * 2.

Numerical methods for solving the resulting systems of differential equations

  * 3.

Distribution of stochastic cost functions

In Section II, we introduce deterministic and stochastic LQ-OCPs and propose
differential equation systems for LQ discretization. For stochastic LQ-OCPs,
we reformulate their cost function and describe the distribution of the
stochastic cost. Section III introduces three numerical methods for solving
proposed differential equation systems. Section IV presents a numerical
experiment comparing the proposed numerical methods, and conclusions are given
in Section V.

##  II Linear-Quadratic Optimal Control Problems

In this section, we introduce deterministic and stochastic LQ-OCPs and
describe the differential equation systems for LQ discretization.

###  II-A Deterministic linear-quadratic optimal control problem

Consider the deterministic LQ-OCP

$\displaystyle\min_{x,u,z,\tilde{z}}\quad$
$\displaystyle\phi=\int_{t_{0}}^{t_{0}+T}l_{c}(\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}\tilde{z}(t)^{\prime}Q_{c}\tilde{z}(t),\end{split}$
(2)

where $Q_{c}=W_{z}^{\prime}W_{z}$ is a semi-positive definite matrix. This
problem is in continuous-time with decision variables $x(t),u(t),z(t)$, and
$\tilde{z}(t)$. The control horizon $T=NT_{s}$ with sampling time $T_{s}$ and
$N\in\mathbb{Z}^{+}$, and $\mathcal{N}=0,1,\ldots,N-1$. We assume piecewise
constant inputs, $u(t)=u_{k}$ and target variables $\bar{z}(t)=\bar{z}_{k}$
for $t_{k}\leq t<t_{k+1}$.

###### Remark 1

Note that the case $\bar{z}(t)=[\bar{x}(t);\,\bar{u}(t)]$,
$C_{c}=\begin{bmatrix}I\\\ 0\end{bmatrix}$, $D_{c}=\begin{bmatrix}0\\\
I\end{bmatrix}$, and $Q_{c}=\begin{bmatrix}Q_{c,xx}&0\\\
0&Q_{c,uu}\end{bmatrix}$ corresponds to

$\begin{split}l_{c}(\tilde{z}(t))&=\frac{1}{2}\left[x(t)-\bar{x}(t)\right]^{\prime}Q_{c,xx}\left[x(t)-\bar{x}(t)\right]\\\
&\qquad+\frac{1}{2}\left[u(t)-\bar{u}(t)\right]^{\prime}Q_{c,uu}\left[u(t)-\bar{u}(t)\right].\end{split}$
(3)

$\blacksquare$

The corresponding discrete-time LQ-OCP is

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

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},$ (5)

where the coefficient in the affine term and the constant term are

$q_{k}=M\bar{z}_{k},\quad\rho_{k}=\int_{t_{k}}^{t_{k+1}}l_{c}(\bar{z}_{k})dt=l_{c}(\bar{z}_{k})T_{s},\quad
k\in\mathcal{N}.$ (6)

######  Proposition 1 (Discretization of the deterministic LQ-OCP)

The system of differential equations

$\displaystyle\dot{A}(t)$ $\displaystyle=A_{c}A(t),\qquad$ $\displaystyle
A(0)=I,$ (7a) $\displaystyle\dot{B}(t)$ $\displaystyle=A(t)B_{c},\qquad$
$\displaystyle B(0)=0,$ (7b) $\displaystyle\dot{Q}(t)$
$\displaystyle=\Gamma(t)^{\prime}Q_{c}\Gamma(t),\qquad$ $\displaystyle
Q(0)=0,$ (7c) $\displaystyle\dot{M}(t)$
$\displaystyle=-\Gamma(t)^{\prime}Q_{c},\qquad$ $\displaystyle M(0)=0,$ (7d)
where
$\Gamma(t)=\begin{bmatrix}C_{c}&D_{c}\end{bmatrix}\begin{bmatrix}A(t)&B(t)\\\
0&I\end{bmatrix},$ (7e)

may be used to compute ($A=A(T_{s})$, $B=B(T_{s})$, $Q=Q(T_{s})$,
$M=M(T_{s})$).

$\blacksquare$

###  II-B Certainty equivalent LQ control for a stochastic system

Consider an initial state and an input noise modelled by the following random
variables,

$\displaystyle{\boldsymbol{x}}(t_{0})$ $\displaystyle\sim
N(\hat{x}_{0},P_{0}),\qquad d{\boldsymbol{\omega}}(t)\sim N_{iid}(0,Idt).$ (8)

The stochastic system can be described as continuous-time linear stochastic
differential equations (SDEs) in the form

$\displaystyle d{\boldsymbol{x}}(t)$
$\displaystyle=\left(A_{c}{\boldsymbol{x}}(t)+B_{c}u(t)\right)dt+G_{c}d{\boldsymbol{\omega}}(t),$
(9a) $\displaystyle{\boldsymbol{z}}(t)$
$\displaystyle=C_{c}{\boldsymbol{x}}(t)+D_{c}u(t).$ (9b)

The corresponding discrete-time stochastic system is

$\displaystyle{\boldsymbol{x}}_{k+1}$
$\displaystyle=A{\boldsymbol{x}_{k}}+Bu_{k}+{\boldsymbol{w}}_{k},$ (10a)
$\displaystyle{\boldsymbol{z}}_{k}$
$\displaystyle=C{\boldsymbol{x}_{k}}+Du_{k},$ (10b) where
$\displaystyle{\boldsymbol{x}}_{0}$ $\displaystyle\sim
N(\hat{x}_{0},P_{0}),\qquad{\boldsymbol{w}}_{k}\sim N_{iid}(0,R_{ww}).$ (10c)

######  Proposition 2 (Discretization of the linear SDE)

The system of differential equations

$\displaystyle\dot{A}(t)$ $\displaystyle=A_{c}A(t),\qquad$ $\displaystyle
A(0)=I,$ (11a) $\displaystyle\dot{B}(t)$ $\displaystyle=A(t)B_{c},\qquad$
$\displaystyle B(0)=0,$ (11b) $\displaystyle\dot{R}_{ww}$
$\displaystyle=\Phi(t)\Phi(t)^{\prime},\qquad$ $\displaystyle R_{ww}(0)=0,$
(11c) where $\Phi(t)=A(t)G_{c},$ (11d)

can be used to compute ($A=A(T_{s})$, $B=B(T_{s})$, $R_{ww}=R_{ww}(T_{s})$).

$\blacksquare$

###  II-C Stochastic linear-quadratic optimal control problem

Consider the stochastic LQ-OCP

$\displaystyle\min_{{\boldsymbol{x}},u,{\boldsymbol{z}},{\tilde{\boldsymbol{z}}}}\quad$
$\displaystyle\psi=E\left\\{\phi=\int_{t_{0}}^{t_{0}+T}l_{c}({\tilde{\boldsymbol{z}}}(t))dt\right\\}$
(12a) $\displaystyle s.t.$ $\displaystyle{\boldsymbol{x}}(t_{0})\sim
N(\hat{x}_{0},P_{0}),$ (12b) $\displaystyle d{\boldsymbol{\omega}}(t)\sim
N_{iid}(0,Idt),$ (12c) $\displaystyle u(t)=u_{k},\quad t_{k}\leq
t<t_{k+1},\,k\in\mathcal{N},$ (12d) $\displaystyle
d{\boldsymbol{x}}(t)=(A_{c}{\boldsymbol{x}}(t)+B_{c}u(t))dt+G_{c}d{\boldsymbol{\omega}}(t),$
(12e) $\displaystyle{\boldsymbol{z}}(t)=C_{c}{\boldsymbol{x}}(t)+D_{c}u(t),$
(12f) $\displaystyle{\bar{z}}(t)=\bar{z}_{k},\quad t_{k}\leq
t<t_{k+1},\,k\in\mathcal{N},$ (12g)
$\displaystyle{\tilde{\boldsymbol{z}}}(t)={\boldsymbol{z}}(t)-\bar{z}(t).$
(12h)

The corresponding discrete-time stochastic LQ-OCP is

$\displaystyle\min_{x,u}\quad$
$\displaystyle\psi=E\left\\{\phi=\sum_{k\in\mathcal{N}}l_{k}({\boldsymbol{x}}_{k},u_{k})+l_{s,k}(\boldsymbol{x}_{k},u_{k})\right\\}$
(13a) $\displaystyle s.t.$ $\displaystyle{\boldsymbol{x}}_{0}\sim
N(\hat{x}_{0},P_{0}),$ (13b) $\displaystyle{\boldsymbol{w}_{k}}\sim
N_{iid}(0,R_{ww}),\qquad\qquad\>k\in\mathcal{N},$ (13c)
$\displaystyle{\boldsymbol{x}}_{k+1}=A{\boldsymbol{x}}_{k}+Bu_{k}+{\boldsymbol{w}}_{k},\qquad
k\in\mathcal{N},$ (13d)

where the stage cost function $l_{k}(\boldsymbol{x}_{k},u_{k})$ is

$l_{k}(\boldsymbol{x}_{k},u_{k})=\frac{1}{2}\begin{bmatrix}\boldsymbol{x}_{k}\\\
u_{k}\end{bmatrix}^{\prime}Q\begin{bmatrix}\boldsymbol{x}_{k}\\\
u_{k}\end{bmatrix}+q_{k}^{\prime}\begin{bmatrix}\boldsymbol{x}_{k}\\\
u_{k}\end{bmatrix}+\rho_{k},$ (14a) and the stochastic stage cost function
$l_{s,k}(\boldsymbol{x}_{k},u_{k})$ is
$l_{s,k}(\boldsymbol{x}_{k},u_{k})=\int_{t_{k}}^{t_{k+1}}\frac{1}{2}\boldsymbol{w}(t)^{\prime}Q_{c,ww}\boldsymbol{w}(t)+\boldsymbol{q}_{s,k}^{\prime}\boldsymbol{w}(t)dt.$
(14b)

$Q$, $q_{k}$, and $\rho_{k}$ in (14a) are identical to the deterministic case.
The state variables and system matrices of $l_{s,k}(\boldsymbol{x}_{k},u_{k})$
are

$\displaystyle\boldsymbol{w}(t)=\int_{0}^{t}A(s)G_{c}d\boldsymbol{\omega}(s),\quad$
$\displaystyle Q_{c,ww}=C_{c}^{\prime}Q_{c,xx}C_{c},$ (15a)
$\displaystyle\tilde{\boldsymbol{z}}_{k}=\Gamma(t)\begin{bmatrix}\boldsymbol{x}_{k}\\\
u_{k}\end{bmatrix}-\bar{z}_{k},$
$\displaystyle\boldsymbol{q}_{s,k}=\begin{bmatrix}C_{c}^{\prime}&0\end{bmatrix}Q_{c}\tilde{\boldsymbol{z}}_{k}.$
(15b)

The expectation of the stochastic LQ-OCP is [3]

$\displaystyle\min_{x,u}\>$
$\displaystyle\psi=\sum_{k\in\mathcal{N}}l_{k}(x_{k},u_{k})+\frac{1}{2}\left[\text{tr}\left(Q\bar{P}_{k}\right)+\int_{t_{k}}^{t_{k+1}}\text{tr}\left(Q_{c,ww}P_{w}\right)dt\right]$
(16a) $\displaystyle s.t.$ $\displaystyle x_{0}=\hat{x}_{0},$ (16b)
$\displaystyle x_{k+1}=Ax_{k}+Bu_{k},\quad k\in\mathcal{N},$ (16c)

where

$\begin{bmatrix}\boldsymbol{x}_{k}\\\ u_{k}\end{bmatrix}\sim
N(m_{k},\bar{P}_{k}),\quad m_{k}=\begin{bmatrix}x_{k}\\\
u_{k}\end{bmatrix},\quad\bar{P}_{k}=\begin{bmatrix}P_{k}&0\\\
0&0\end{bmatrix},$ (17a) $P_{k+1}=AP_{k}A^{\prime}+R_{ww},\qquad
P_{w}=\text{Cov}\left(\boldsymbol{w}(t)\right).$ (17b)

######  Proposition 3 (Discretization of the stochastic LQ-OCP)

The system of differential equations

$\displaystyle\dot{A}(t)$ $\displaystyle=A_{c}A(t),\qquad$ $\displaystyle
A(0)=I,$ (18a) $\displaystyle\dot{B}(t)$ $\displaystyle=A(t)B_{c},\qquad$
$\displaystyle B(0)=0,$ (18b) $\displaystyle\dot{R}_{ww}(t)$
$\displaystyle=\Phi(t)\Phi(t)^{\prime},\qquad$ $\displaystyle R_{ww}(0)=0,$
(18c) $\displaystyle\dot{Q}(t)$
$\displaystyle=\Gamma(t)^{\prime}Q_{c}\Gamma(t),$ $\displaystyle Q(0)=0,$
(18d) $\displaystyle\dot{M}(t)$
$\displaystyle=-\Gamma(t)^{\prime}Q_{c},\qquad$ $\displaystyle M(0)=0,$ (18e)
where $\displaystyle\Phi(t)$
$\displaystyle=A(t)G_{c},\qquad\Gamma(t)=\begin{bmatrix}C_{c}&D_{c}\end{bmatrix}\begin{bmatrix}A(t)&B(t)\\\
0&I\end{bmatrix},$ (18f)

can be used to compute ($A=A(T_{s})$, $B=B(T_{s})$, $R_{ww}=R_{ww}(T_{s})$,
$Q=Q(T_{s})$, $M=M(T_{s})$).

$\blacksquare$

######  Proposition 4 (Distribution of the stochastic costs)

Using Euler-Maruyama (EM) discretization, we can reformulate (13a) into
isolated stochastic form as

$\begin{split}\phi&=\sum_{k\in\mathcal{N}}l_{k}(\boldsymbol{x}_{k},u_{k})+l_{s,k}(\boldsymbol{x}_{k},u_{k})\\\
&=\frac{1}{2}\begin{bmatrix}\boldsymbol{x}_{0}\\\
\boldsymbol{W}_{\mathcal{N}}\end{bmatrix}^{\prime}Q_{\mathcal{N}}\begin{bmatrix}\boldsymbol{x}_{0}\\\
\boldsymbol{W}_{\mathcal{N}}\end{bmatrix}+q_{\mathcal{N}}^{\prime}\begin{bmatrix}\boldsymbol{x}_{0}\\\
\boldsymbol{W}_{\mathcal{N}}\end{bmatrix}+\rho_{\mathcal{N}},\end{split}$
(19a) where $\begin{bmatrix}\boldsymbol{x}_{0}\\\
\boldsymbol{W}_{\mathcal{N}}\end{bmatrix}\sim
N(\bar{m},\bar{P}),\quad\bar{m}=\begin{bmatrix}x_{0}\\\
0\end{bmatrix},\quad\bar{P}=\begin{bmatrix}P_{0}&0\\\ 0&P_{w}\end{bmatrix}.$
(19b)

$\boldsymbol{W}_{\mathcal{N}}$ is a vector of sub-sampling random variables
over the horizon, and its covariance is $P_{w}=\text{diag}([I\delta t,I\delta
t,\ldots,I\delta t])$.

Based on [9], the reformulated stochastic cost follows a generalized
$\chi^{2}$ distribution, and its expectation and variance can be expressed as

$\displaystyle
E\left\\{\phi\right\\}=\frac{1}{2}\bar{m}^{\prime}Q_{\mathcal{N}}\bar{m}+q_{\mathcal{N}}^{\prime}\bar{m}+\rho_{\mathcal{N}}+\frac{1}{2}\text{tr}\left(Q_{\mathcal{N}}\bar{P}\right),$
(20a)
$\displaystyle\begin{split}&V\left\\{\phi\right\\}=q_{\mathcal{N}}^{\prime}\bar{P}q_{\mathcal{N}}+2\bar{m}^{\prime}Q_{\mathcal{N}}\bar{P}q_{\mathcal{N}}+\bar{m}^{\prime}Q_{\mathcal{N}}\bar{P}Q_{\mathcal{N}}\bar{m}\\\
&\qquad\qquad+\frac{1}{2}\text{tr}(Q_{\mathcal{N}}\bar{P}Q_{\mathcal{N}}\bar{P}).\end{split}$
(20b)

The stochastic cost $\phi$ in (12a) follows a generalized $\chi^{2}$
distribution with mean and variance given by (20) when taking the limit of
integration steps $n\rightarrow\infty$. See Appendix Distribution of
stochastic costs for more details.

$\blacksquare$

##  III Numerical Methods of LQ Discretization

This section introduces numerical methods for the discretization of
continuous-time LQ-OCPs.

###  III-A Ordinary differential equation methods

Table I: $\Lambda$ and $\Theta$ of the ODE method with discretization methods Methods |  $\Lambda$ |  $\Theta$
---|---|---
Expl. Euler |  $I+hA_{c}$ |  $I$
Impl. Euler |  $(I-hA_{c})^{-1}$ |  $(I-hA_{c})^{-1}$
Expl. Trape. |  $I+hA_{c}+0.5h^{2}A_{c}^{2}$ |  $I+0.5hA_{c}$
Impl. Trape. |  $(I+0.5hA_{c})^{-1}$ $(I-0.5hA_{c})$ |  $(I-0.5hA_{c})^{-1}$
ESDIRK34 |  $(I-0.44hA_{c})^{-3}$ $(I-0.31hA_{c}-$ $0.24h^{2}A_{c}^{2})$ |  $(I-0.44hAc)^{-3}$ $(I-0.81hAc+$ $0.08h^{2}Ac^{2})$
Classic RK4 |  $I+hA_{c}+0.5h^{2}A_{c}^{2}+$ $0.17h^{3}A_{c}^{3}+$ $0.04h^{4}A_{c}^{4}$ |  $I+0.5hA_{c}+0.17h^{2}A_{c}^{2}+$ $0.04h^{3}A_{c}^{3}$

Consider an s-stage ODE method with $N\in\mathbb{Z}^{+}$ integration steps and
the time step $h=\frac{T_{s}}{N}$. We compute ($A$, $B$, $R_{ww}$, $Q$, $M$)
as

$\displaystyle A_{k+1}$ $\displaystyle=\Lambda A_{k},\>$ $\displaystyle
k\in\mathcal{N},$ (21a) $\displaystyle B_{k+1}$ $\displaystyle=B_{k}+\Theta
A_{k}\bar{B}_{c},\>$ $\displaystyle k\in\mathcal{N},$ (21b)
$\displaystyle\Gamma_{k+1}$ $\displaystyle=\Omega\Gamma_{k},\>$ $\displaystyle
k\in\mathcal{N},$ (21c) $\displaystyle M_{k+1}$
$\displaystyle=M_{k}+\Gamma_{k}^{\prime}\bar{M}_{c},\>$ $\displaystyle
k\in\mathcal{N},$ (21d) $\displaystyle Q_{k+1}$
$\displaystyle=Q_{k}+\Gamma_{k}^{\prime}\bar{Q}_{c}\Gamma_{k},\>$
$\displaystyle k\in\mathcal{N},$ (21e) $\displaystyle R_{ww,k+1}$
$\displaystyle=R_{ww,k}+\sum_{i=1}^{s}b_{i}\Lambda_{i}A_{k}\bar{R}_{ww,c}A_{k}^{\prime}\Lambda_{i}^{\prime},\>$
$\displaystyle k\in\mathcal{N},$ (21f)

where $\bar{B}_{c}=hB_{c}$, $H=\begin{bmatrix}C_{c}&D_{c}\end{bmatrix}$,
$\bar{M}_{c}=-h\sum_{i=1}^{s}b_{i}\Omega_{i}^{\prime}H^{\prime}Q_{c}$,
$\bar{Q}_{c}=h\sum_{i=1}^{s}b_{i}\Omega_{i}^{\prime}H^{\prime}Q_{c}H\Omega_{i}$,
and $\bar{R}_{ww,c}=hG_{c}G_{c}^{\prime}$ are constant matrices. $\Lambda$,
$\Theta$, $\Omega$ are functions of the coefficients, $a_{i,j}$ and $b_{i}$
for $i=1,2,\ldots,s$ and $j=1,2,\ldots,s$ in the Butcher tableau of the ODE
method. They are computed as

$\displaystyle A_{k,i}$
$\displaystyle=A_{k}+h\sum_{j=1}^{s}a_{i,j}\dot{A}_{k,j}=\Lambda_{i}A_{k},$
(22a) $\displaystyle B_{k,i}$
$\displaystyle=B_{k}+h\sum_{j=1}^{s}a_{i,j}\dot{B}_{k,j}=B_{k}+\Theta_{i}A_{k}\bar{B}_{c},$
(22b) $\displaystyle\Gamma_{k,i}$
$\displaystyle=\Omega_{i}\Gamma_{k}=\begin{bmatrix}\Lambda_{i}&\Theta_{i}\bar{B}_{c}\\\
0&I\end{bmatrix}\Gamma_{k},$ (22c) $\displaystyle\Lambda$
$\displaystyle=I+h\sum_{i=1}^{s}b_{i}A_{c}\Lambda_{i},$ (22d)
$\displaystyle\Theta$ $\displaystyle=\sum_{i=1}^{s}b_{i}\Lambda_{i},$ (22e)
$\displaystyle\Omega$
$\displaystyle=\begin{bmatrix}\Lambda&\Theta\bar{B}_{c}\\\ 0&I\end{bmatrix},$
(22f)

where $\Lambda_{i}$, $\Theta_{i}$, $\Omega_{i}$ are coefficients of stage
variables $A_{k,i}$, $B_{k,i}$, $\Gamma_{k,i}$.

Consequently, the differential equation systems $A(T_{s})=A_{N}$,
$B(T_{s})=B_{N}$, $R_{ww}(T_{s})=R_{ww,N}$, $Q(T_{s})=Q_{N}$,
$M(T_{s})=M_{N}$, has constant coefficients $\Lambda_{i}$, $\Theta_{i}$,
$\Omega_{i}$, $\Lambda$, $\Theta$ and $\Omega$ when using fixed-time-step ODE
methods. The constant coefficients can be computed offline. Table I describes
$\Lambda$ and $\Theta$ for different discretization methods. Algorithm 1
presents the ODE methods for LQ discretization.

Algorithm 1 ODE method for LQ Discretization

Input: $(A_{c},B_{c},G_{c},C_{c},D_{c},Q_{c},T_{s},h)$

Output: $(A(T_{s}),B(T_{s}),C,D,Q(T_{s}),M(T_{s}),R_{ww}(T_{s}))$

Set initial states

($k=0$, $A_{k}=I$, $B_{k}=0$, $Q_{k}=0$, $M_{k}=0$, $R_{ww,k}=0$)

Use (22) to compute ($\Lambda_{i}$, $\Theta_{i}$, $\Omega_{i}$, $\Lambda$,
$\Theta$, $\Omega$)

Compute integration steps $N=\frac{T_{s}}{h}$

while $k<N$ do

Use (21) to update $(A_{k},B_{k},\Gamma_{k},Q_{k},M_{k},R_{ww,k})$

Set $k\leftarrow k+1$

end while

Get system matrices
$(A(T_{s})=A_{k},B(T_{s})=B_{k},C=C_{c},D=D_{c},Q(T_{s})=Q_{k},M(T_{s})=M_{k},R_{ww}(T_{s})=R_{ww,k})$

###  III-B Matrix exponential method

The matrix exponential method describes the LQ discretization by three matrix
exponential problems

$\begin{split}\begin{bmatrix}\Phi_{1,11}&\Phi_{1,12}\\\
0&\Phi_{1,22}\end{bmatrix}&=\text{exp}\left(\begin{bmatrix}-H^{\prime}&\bar{Q}_{c}\\\
0&H\end{bmatrix}t\right),\\\ \end{split}$ (23a)
$\begin{split}\begin{bmatrix}I&\Phi_{2,12}\\\
0&\Phi_{2,22}\end{bmatrix}&=\text{exp}\left(\begin{bmatrix}0&I\\\
0&H^{\prime}\end{bmatrix}t\right),\\\ \end{split}$ (23b)
$\begin{split}\begin{bmatrix}\Phi_{3,11}&\Phi_{3,12}\\\
0&\Phi_{3,22}\end{bmatrix}&=\text{exp}\left(\begin{bmatrix}-A_{c}&\bar{G}_{c}\\\
0&A_{c}^{\prime}\end{bmatrix}t\right),\\\ \end{split}$ (23c)

where

$\displaystyle H$ $\displaystyle=\begin{bmatrix}A_{c}&B_{c}\\\
0&0\end{bmatrix},\qquad$
$\displaystyle\bar{M}_{c}=-\begin{bmatrix}C_{c}&D_{c}\end{bmatrix}^{\prime}Q_{c},$
(24a) $\displaystyle\bar{Q}_{c}$
$\displaystyle=-\bar{M}_{c}\begin{bmatrix}C_{c}&D_{c}\end{bmatrix},\qquad$
$\displaystyle\bar{G}_{c}=G_{c}G_{c}^{\prime}.$ (24b)

The elements of matrix exponential problems are

$\displaystyle\begin{split}\Phi_{1,22}&=\Gamma(t)=\begin{bmatrix}A(t)&B(t)\\\
0&I\end{bmatrix},\end{split}$ (25a)
$\displaystyle\begin{split}\Phi_{1,12}&=\Gamma(-t)^{\prime}\int_{0}^{t}\Gamma(\tau)^{\prime}\bar{Q}_{c}\Gamma(\tau)d\tau,\end{split}$
(25b)
$\displaystyle\begin{split}\Phi_{2,12}&=\int_{0}^{t}\Gamma(\tau)^{\prime}d\tau,\end{split}$
(25c) $\displaystyle\begin{split}\Phi_{3,22}&=A(t)^{\prime},\end{split}$ (25d)
$\displaystyle\begin{split}\Phi_{3,12}&=A(-t)^{\prime}\int^{t}_{0}A(\tau)\bar{G}_{c}A(\tau)^{\prime}d\tau.\end{split}$
(25e)

Set $t=T_{s}$, we can compute differential equations $A$, $B$, $R_{ww}$, $Q$
and $M$ as

$\displaystyle A(T_{s})$ $\displaystyle=\Phi_{1,22}(1:n_{x},1:n_{x}),$ (26a)
$\displaystyle B(T_{s})$
$\displaystyle=\Phi_{1,22}(1:n_{x},n_{x}+1:\text{end}),$ (26b) $\displaystyle
Q(T_{s})$ $\displaystyle=\Phi_{1,22}^{\prime}\Phi_{1,12},$ (26c)
$\displaystyle M(T_{s})$ $\displaystyle=\Phi_{2,12}\bar{M}_{c},$ (26d)
$\displaystyle R_{ww}(T_{s})$ $\displaystyle=\Phi_{3,22}^{\prime}\Phi_{3,12}.$
(26e)

The matrix exponential method is inspired by formulas from [28, 29, 34].

###  III-C Step-doubling method

Consider an s-stage ODE method with $N\in\mathbb{Z}^{+}$ integration steps and
the time step $h=\frac{T_{s}}{N}$. The matrices

$\displaystyle\tilde{A}(N)$ $\displaystyle=\Lambda^{N},$
$\displaystyle\tilde{A}(1)=\Lambda,$ (27a) $\displaystyle\tilde{B}(N)$
$\displaystyle=\displaystyle\sum_{i=0}^{N-1}\Lambda^{i},$
$\displaystyle\tilde{B}(1)=I,$ (27b) $\displaystyle\tilde{\Gamma}(N)$
$\displaystyle=\Omega^{N},$ $\displaystyle\tilde{\Gamma}(1)=I_{xu},$ (27c)
$\displaystyle\tilde{M}(N)$
$\displaystyle=\displaystyle\sum_{i=0}^{N-1}\left(\Gamma^{i}\right)^{\prime},$
$\displaystyle\tilde{M}(1)=I_{xu},$ (27d) $\displaystyle\tilde{Q}(N)$
$\displaystyle=\displaystyle\sum_{i=0}^{N-1}\left(\Gamma^{i}\right)^{\prime}\bar{Q}_{c}\left(\Gamma^{i}\right),$
$\displaystyle\tilde{Q}(1)=\bar{Q}_{c},$ (27e) $\displaystyle\tilde{R}(N)$
$\displaystyle=\displaystyle\sum_{i=0}^{N-1}\left(A^{i}\right)\bar{R}_{ww,c}\left(A^{i}\right)^{\prime},$
$\displaystyle\tilde{R}(1)=\bar{R}_{ww,c},$ (27f)

can be used to solve $(A,B,M,Q,R_{ww})$

$\displaystyle A(T_{s})$ $\displaystyle=\tilde{A}(N),$ (28a) $\displaystyle
B(T_{s})$ $\displaystyle=\Theta\tilde{B}(N)\bar{B}_{c},$ (28b) $\displaystyle
M(T_{s})$ $\displaystyle=\tilde{M}(N)\bar{M}_{c},$ (28c) $\displaystyle
Q(T_{s})$ $\displaystyle=\tilde{Q}(N),$ (28d) $\displaystyle R_{ww}(T_{s})$
$\displaystyle=\sum_{i=1}^{s}b_{i}\Lambda_{i}\tilde{R}(N)\Lambda_{i}^{\prime},$
(28e)

where $\Lambda$, $\Theta$, $\Omega$, $\bar{B}_{c}$, $\bar{M}_{c}$ and
$\bar{Q}_{c}$ are the same as in the ODE method case.
$I_{xu}\in\mathbb{R}^{n_{xu}\times n_{xu}}$ is an identity matrix with the
size $n_{xu}=n_{x}+n_{u}$.

Figure 1: The exponential of a matrix $A_{c}$. The blue line is the true
result computed with expm() in MATLAB. The red dots are the results of the
step-doubling method.

Fig. 1 shows an example of $\tilde{A}(t)=e^{tA_{c}}$, where the step-doubling
method (red dots) uses the $\frac{n}{2}^{th}$ step result as the initial state
to compute the double step’s result $\tilde{A}(n)$. We then get the step-
doubling expression as

$\tilde{A}(1)\rightarrow\tilde{A}(2)\rightarrow\tilde{A}(4)\rightarrow\ldots\rightarrow\tilde{A}(\frac{N}{4})\rightarrow\tilde{A}(\frac{N}{2})\rightarrow\tilde{A}(N),$
(29a) where
$\tilde{A}(n)=\tilde{A}(\frac{n}{2})\tilde{A}(\frac{n}{2}),\qquad\qquad
n\in[2,4,\ldots,\frac{N}{2},N].$ (29b)

Eq. (29b) is inspired by the scaling and squaring algorithm for solving matrix
exponential problem presented in [1, 2, 17]. We apply the same idea for other
differential equations. Table II describes the step-doubling expressions for
($\tilde{A}$, $\tilde{B}$, $\tilde{\Gamma}$, $\tilde{Q}$, $\tilde{M}$,
$\tilde{R}_{ww}$). The step-doubling method takes only $j$ steps to get the
same result as the ODE method with $N=2^{j}$ integration steps. Algorithm 2
describes the step-doubling method for LQ discretization.

Table II: Numerical expressions of the step-doubling method ODEs |  ODE expression |  Step-doubling expression
---|---|---
$\tilde{A}(N)$ |  $\Lambda^{N}$ |  $\tilde{A}(\frac{N}{2})\tilde{A}(\frac{N}{2})$
$\tilde{B}(N)$ |  $\displaystyle\sum_{i=0}^{N-1}\Lambda^{i}$ |  $\tilde{B}(\frac{N}{2})\left(I+\tilde{A}(\frac{N}{2})\right)$
$\tilde{\Gamma}(N)$ |  $\Omega^{N}$ |  $\tilde{\Gamma}(\frac{N}{2})\tilde{\Gamma}(\frac{N}{2})$
$\tilde{M}(N)$ |  $\displaystyle\sum_{i=0}^{N-1}\Omega^{i}\bar{M}_{c}$ |  $\tilde{M}(\frac{N}{2})\left(I+\tilde{\Gamma}(\frac{N}{2})^{\prime}\right)$
$\tilde{Q}(N)$ |  $\displaystyle\sum_{i=0}^{N-1}\left(\Gamma^{i}\right)^{\prime}\bar{Q}_{c}\left(\Gamma^{i}\right)$ |  $\tilde{Q}(\frac{N}{2})+\tilde{\Gamma}(\frac{N}{2})^{\prime}\tilde{Q}(\frac{N}{2})\tilde{\Gamma}(\frac{N}{2})$
$\tilde{R}(N)$ |  $\displaystyle\sum_{i=0}^{N-1}\left(A^{i}\right)\bar{R}_{ww,c}\left(A^{i}\right)^{\prime}$ |  $\tilde{R}(\frac{N}{2})+\tilde{A}(\frac{N}{2})\tilde{R}(\frac{N}{2})\tilde{A}(\frac{N}{2})^{\prime}$
Algorithm 2 Step-doubling method for LQ Discretization

Input: $(A_{c},B_{c},G_{c},C_{c},D_{c},Q_{c},T_{s},j)$

Output: $(A(T_{s}),B(T_{s}),C,D,Q(T_{s}),M(T_{s}),R_{ww}(T_{s}))$

Compute the number of step $N=2^{j}$

Compute the step size $h=\frac{T_{s}}{N}$

Use (22) to compute ($\Lambda_{i}$, $\Theta_{i}$, $\Omega_{i}$, $\Lambda$,
$\Theta$, $\Omega$)

Set initial states of step doubling matrices ($i=1,\tilde{A}(i)=\Lambda$,
$\tilde{B}(i)=I$, $\tilde{Q}(i)=\bar{Q}_{c}$, $\tilde{M}(i)=I_{xu}$,
$\tilde{R}(i)=\bar{R}_{ww,c}$)

while $i\leq j$ do

Use equations from Table II to update ($\tilde{\Gamma}(i)$,$\tilde{M}(i)$,
$\tilde{Q}(i)$,$\tilde{R}(i)$)

Use equations from Table II to update ($\tilde{A}(i)$, $\tilde{B}(i)$)

Set $i=i+1$

end while

Use (28) to compute $(A(T_{s}),B(T_{s}),Q(T_{s}),M(T_{s}),R_{ww}(T_{s}))$

##  IV Numerical Experiments

Consider a continuous-time stochastic LQ-OCP with system matrices

$A_{c}=\begin{bmatrix}-49&24\\\ -64&31\end{bmatrix},\quad
B_{c}=\begin{bmatrix}2&0.5\\\ 1&3\end{bmatrix},\quad
G_{c}=\begin{bmatrix}0.1&0.0\\\ 0.0&0.1\end{bmatrix}.$ (30)

The system output matrices are $C_{c}=[1.0,1.0]$, $D_{c}=[0.0,0.0]$, and the
system references and system inputs are

$\displaystyle\bar{z}(t)=\bar{z}_{k}=3.0,$ $\displaystyle t_{k}\leq t\leq
t_{k+1},\quad k\in\mathcal{N},$ (31a)
$\displaystyle\bar{u}(t)=\bar{u}_{k}=0.0,$ $\displaystyle t_{k}\leq t\leq
t_{k+1},\quad k\in\mathcal{N},$ (31b) $\displaystyle
u(t)=u_{k}=\begin{bmatrix}1.0&1.0\end{bmatrix}^{T},$ $\displaystyle t_{k}\leq
t\leq t_{k+1},\quad k\in\mathcal{N}.$ (31c)

The weights are $Q_{c,xx}=1.0$ and $Q_{c,uu}=\text{diag}([1.0,1.0])$, the
sampling time $T_{s}=1.0$, and initial state vector is
$\boldsymbol{x}_{0}=\begin{bmatrix}0.0&1.0\end{bmatrix}^{T}$ with the
covariance $P_{0}=\text{diag}([0.1,0.1])$.

In this section, we test and compare the proposed numerical methods and
investigate the distribution of the stochastic costs via Monte Carlo
simulations.

###  IV-A Discretization of LQ-OCP

Figure 2: The error and CPU time of the ODE methods and the step-doubling
methods with different discretization methods. The error is
$e(i)=|i(T_{s})-i(N)|$ for $i\in[A,B,R_{ww},M,Q]$, where $i(T_{s})$ is the
true result from the matrix exponential method.

Fig. 2 describes the error and CPU time of ODE and step-doubling methods. The
true solution of ($A,B,R_{ww},Q,M$) is calculated using the matrix exponential
method. The results of the step-doubling method (dot plots) have the same
error as the results of the ODE method (line plots). All methods have the
correct convergence order (indicated by dashed lines). In bar plots, the CPU
time of the ODE method increases as the integration steps and the stages of
the discretization method increase. However, the CPU time of the step-doubling
method is stable at around 0.6 ms.

Table III: CPU time and error of the scenario using classic RK4 with $N=2^{8}$ | Unit | Matrix Exp. | ODE Method | Step-doubling
---|---|---|---|---
$e(A)$ | [-] | - | $7.49\cdot 10^{-12}$ | $7.49\cdot 10^{-12}$
$e(B)$ | [-] | - | $8.33\cdot 10^{-12}$ | $8.33\cdot 10^{-12}$
$e(R_{ww})$ | [-] | - | $9.73\cdot 10^{-11}$ | $9.73\cdot 10^{-11}$
$e(M)$ | [-] | - | $1.25\cdot 10^{-11}$ | $1.25\cdot 10^{-11}$
$e(Q)$ | [-] | - | $2.03\cdot 10^{-13}$ | $2.03\cdot 10^{-13}$
CPU Time | [ms] | 0.74 | 9.5 | 0.68

Table III describes the error and CPU time of ODE and step-doubling methods
with the classic RK4 method applied with the integration step $N=2^{8}$. The
step-doubling method has the same error as the ODE method. The ODE method is
the slowest among the three methods and takes 9.5 ms, while the matrix
exponential and the step-doubling methods spend 0.74 ms and 0.68 ms,
respectively.

###  IV-B Distribution of stochastic LQ-OCP

Figure 3: The likelihood of the cost functions of continuous-time and
discrete-time stochastic LQ-OCPs with 30000 Monte Carlo simulations. The
Cont., Disc. EM. indicate the continuous-time, discrete-time, and EM
reformulated stochastic costs, respectively. $E\\{\phi\\}$-Analy. is the
analytic expectation described by (16), where the continuous-time element
$\text{tr}(Q_{c,ww}P_{w})$ is solved using the EM method with $N=2^{8}$.

Fig. 3 shows the likelihood of the continuous-time stochastic costs (12a),
discrete stochastic costs (13a), and EM reformulated stochastic costs (19a)
via 30,000 Monte Carlo simulations. We apply the EM method to solve these
stochastic cost functions with $N=2^{8}$ integration steps.

The simulation results show that the continuous-time cost has the same
distribution as the discrete and EM reformulated stochastic costs. There is an
offset between the analytic expectation ($E(\phi)=6.36$, obtained by (16)) and
the other numerical expectations (6.47 for continuous- and discrete-time cases
and 6.46 for the EM case). We consider it reasonable to have numerical errors
since we cannot take the limit of $N\rightarrow\infty$ in experiments.

##  V Conclusions

In this paper, we have discussed the discretization of both deterministic and
stochastic LQ-OCPs and proposed three numerical methods. In our propositions,
LQ discretization is converted into explicit, neat differential equation
systems. We further extend the problem from the deterministic to the
stochastic case and illustrate the stochastic cost adheres to a generalized
$\chi^{2}$ distribution. The proposed numerical methods are tested and
compared in the numerical experiment, and its results indicate: 1) the step-
doubling method is the fastest among the three methods while retaining the
same accuracy and convergence order as the ODE method, 2) the discrete-time
LQ-OCP derived by the proposed numerical methods is equivalent to the original
problem in both the deterministic and stochastic cases.

##  VI APPENDIX

## Distribution of stochastic costs

To evaluate the distribution of the costs, consider the EM discretization of
the stochastic system with a fine time step $\delta t=\frac{T_{s}}{n}$,
$t_{i}=i\delta t$ for $i=1,\ldots,n$ is

$x_{k,i}=\overbrace{(I+\delta
tA_{c})^{i}}^{A_{i}}x_{k}+\overbrace{\left(\sum_{j=0}^{i-1}A(j)\delta
tB_{c}\right)}^{B_{i}}u_{k}+G_{i}\boldsymbol{w}_{k},$ (32a) where
$\displaystyle G_{i}$
$\displaystyle=\begin{bmatrix}A_{i-1}G_{c}&A_{i-2}G_{c}&\ldots&0_{n-i}\end{bmatrix},$
(32b) $\displaystyle\boldsymbol{w}_{k}$
$\displaystyle=\begin{bmatrix}\Delta\boldsymbol{w}_{k,1}^{\prime}&\Delta\boldsymbol{w}_{k,2}^{\prime}&\ldots&\Delta\boldsymbol{w}_{k,n}^{\prime}\end{bmatrix}^{\prime}.$
(32c)

The EM expression for the extended state vector $[x_{k};u_{k}]$ is

$\begin{bmatrix}\boldsymbol{x}_{k}\\\
u_{k}\end{bmatrix}=\begin{bmatrix}A_{n}^{k}\\\
0\end{bmatrix}\boldsymbol{x}_{0}+\begin{bmatrix}\Theta_{u,k}\\\
I_{u,k}\end{bmatrix}U_{\mathcal{N}}+\begin{bmatrix}\Theta_{w,k}\\\
0\end{bmatrix}\boldsymbol{W}_{\mathcal{N}},$ (33a) where
$\boldsymbol{W}_{\mathcal{N}}$ and $U_{\mathcal{N}}$ are vectors of the random
$\boldsymbol{w}_{k}=I_{w,k}\boldsymbol{W}_{\mathcal{N}}$ and the input
$u_{k}=I_{u,k}U_{\mathcal{N}}$ over the horizon $\mathcal{N}$, and
$\displaystyle\Theta_{u,k}$
$\displaystyle=\begin{bmatrix}A^{k-1}_{n}B_{n}&A^{k-2}_{n}B_{n}&\ldots&0_{\mathcal{N}-k}\end{bmatrix},$
(33b) $\displaystyle\Theta_{w,k}$
$\displaystyle=\begin{bmatrix}A^{k-1}_{n}G_{n}&A^{k-2}_{n}G_{n}&\ldots&0_{\mathcal{N}-k}\end{bmatrix}.$
(33c)

The corresponding EM expression of the discrete stochastic cost function
$\phi$ is

$\begin{split}\phi&=\sum_{k\in\mathcal{N}}l_{k}(\boldsymbol{x}_{k},u_{k})+l_{s,k}(\boldsymbol{x}_{k},u_{k})\\\
&=\sum_{k\in\mathcal{N}}\frac{1}{2}\begin{bmatrix}\boldsymbol{x}_{k}\\\
u_{k}\end{bmatrix}^{\prime}Q\begin{bmatrix}\boldsymbol{x}_{k}\\\
u_{k}\end{bmatrix}+q_{k}^{\prime}\begin{bmatrix}\boldsymbol{x}_{k}\\\
u_{k}\end{bmatrix}+\rho_{k},\end{split}$ (34)

where

$\displaystyle q_{k}$ $\displaystyle=M\bar{z}_{k}+\Omega\boldsymbol{w}_{k},$
(35a)
$\displaystyle\begin{split}\rho_{k}&=\sum_{i=1}^{n}\frac{1}{2}\boldsymbol{w}_{k}^{\prime}(G_{i}^{\prime}\bar{Q}_{c,ww}G_{i})\boldsymbol{w}_{k}+(-\delta
tG_{i}^{\prime}\bar{Q}_{c}\bar{z}_{k})^{\prime}\boldsymbol{w}_{k}\\\
&\qquad+\frac{1}{2}\bar{z}_{k}^{\prime}\bar{Q}_{c}\bar{z}_{k}.\end{split}$
(35b) $\displaystyle\Gamma_{i}$
$\displaystyle=\begin{bmatrix}C_{c}&D_{c}\end{bmatrix}\begin{bmatrix}A_{i}&B_{i}\\\
0&I\end{bmatrix},$ (35c) $\displaystyle\Omega$
$\displaystyle=\sum_{i=1}^{n}\delta
t\Gamma_{i}^{\prime}Q_{c}\begin{bmatrix}C_{c}^{\prime}&0\end{bmatrix}^{\prime}G_{i},$
(35d)

Using (33a) to substitute $[x_{k};u_{k}]$, we reformulate the quadratic
problem into isolated stochastic form (36a) with deterministic weights
$Q_{\mathcal{N}}$, $q_{\mathcal{N}}$ and $\rho_{\mathcal{N}}$.

$\displaystyle\phi=\frac{1}{2}\begin{bmatrix}\boldsymbol{x}_{0}\\\
\boldsymbol{W}_{\mathcal{N}}\end{bmatrix}^{\prime}Q_{\mathcal{N}}\begin{bmatrix}\boldsymbol{x}_{0}\\\
\boldsymbol{W}_{\mathcal{N}}\end{bmatrix}+q_{\mathcal{N}}^{\prime}\begin{bmatrix}\boldsymbol{x}_{0}\\\
\boldsymbol{W}_{\mathcal{N}}\end{bmatrix}+\rho_{\mathcal{N}},$ (36a)
$\displaystyle\begin{bmatrix}\boldsymbol{x}_{0}\\\
\boldsymbol{W}_{\mathcal{N}}\end{bmatrix}\sim
N(\bar{m},\bar{P}),\quad\bar{m}=\begin{bmatrix}x_{0}\\\
0\end{bmatrix},\quad\bar{P}=\begin{bmatrix}P_{0}&0\\\ 0&P_{w}\end{bmatrix},$
(36b)

where the state vector of the reformulated cost function is normally
distributed, and $\boldsymbol{W}_{\mathcal{N}}$ has the covariance
$P_{w}=\text{diag}(I\delta t,I\delta t,\ldots,I\delta t)$.

Based on the theory of integrating the normal in the quadratic domain [9], the
stochastic cost $\phi$ follows a generalized $\chi^{2}$ distribution. The
quadratic-form of its expectation and variance can be computed as (20a) and
(20b).

The original cost (12a) is equivalent to the discrete cost when taking the
limit

$\displaystyle\int^{t_{0}+T}_{t_{0}}l_{c}(\tilde{z}(t))dt=\lim_{n\rightarrow\infty}\sum^{n}_{j=1}l_{c}(\tilde{z}(t_{0}+j\delta
t))\delta t,\quad\delta t=\frac{T}{n}.$ (37)

Thus, the stochastic cost $\phi$ of the continuous-time LQ-OCP is a
generalized $\chi^{2}$-distribution variable, with mean and variance taken in
the limit of (20).

## References

  * [1] A. H. Al-Mohy and N. J. Higham, “A new scaling and squaring algorithm for the matrix exponential,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 31, no. 3, pp. 970–989, 2010.
  * [2] ——, “Computing the action of the matrix exponential, with an application to exponential integrators,” _SIAM Journal on Scientific Computing_ , vol. 33, no. 2, pp. 488–511, 2011.
  * [3] K. J. Åström, _Introduction to stochastic control theory_. Academic Press New York, 1970.
  * [4] K. J. Åström and B. Wittenmark, _Computer-Controlled Systems: Theory and Design, Third Edition_ , ser. Dover Books on Electrical Engineering. Dover Publications, 2011.
  * [5] A. E. Bryson and Y. C. Ho, _Applied Optimal Control: Optimization, Estimation, and Control (1st ed.)_. New York, USA: Routledge, 1975.
  * [6] A. Capolei, B. Foss, and J. B. Jørgensen, “Profit and risk measures in oil production optimization,” _IFAC-PapersOnLine_ , vol. 48, no. 6, pp. 214–220, 2015, 2nd IFAC Workshop on Automatic Control in Offshore Oil and Gas Production OOGP 2015.
  * [7] B. A. Conway, _Spacecraft Trajectory Optimization_ , ser. Cambridge Aerospace Series. Cambridge University Press, 2010.
  * [8] J. Cullum, “Discrete approximations to continuous optimal control problems,” _SIAM Journal on Control_ , vol. 7, no. 1, pp. 32–49, 1969.
  * [9] A. Das and W. S. Geisler, “A method to integrate and classify normal distributions,” _Journal of Vision_ , vol. 21, no. 10, pp. 1–1, 09 2021\.
  * [10] A. L. Dontchev and W. W. Hager, “The euler approximation in state constrained optimal control,” _Mathematics of Computation_ , vol. 70, no. 233, pp. 173–203, 2001.
  * [11] 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.
  * [12] 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.
  * [13] R. Goebel and M. Subbotin, “Continuous time linear quadratic regulator with control constraints via convex duality,” _IEEE Transactions on Automatic Control_ , vol. 52, no. 5, pp. 886–892, 2007.
  * [14] G. C. Goodwin, J. C. Aguero, M. E. Cea Garridos, M. E. Salgado, and J. I. Yuz, “Sampling and sampled-data models: The interface between the continuous world and digital algorithms,” _IEEE Control Systems Magazine_ , vol. 33, no. 5, pp. 34–53, 2013.
  * [15] L. Han, M. Camlibel, J. Pang, and W. Heemels, “Convergence of discrete-time approximations of constrained linear-quadratic optimal control problems,” in _49th IEEE Conference on Decision and Control (CDC)_ , 2010, pp. 5210–5215.
  * [16] E. Hendricks, P. H. Sørensen, and O. Jannerup, _Linear Systems Control: Deterministic and Stochastic Methods_. Berlin, German: Springer, 2008.
  * [17] N. J. Higham, “The scaling and squaring method for the matrix exponential revisited,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 26, no. 4, pp. 1179–1193, 2005.
  * [18] S. Hørsholt, “Operational optimization of oil reservoirs in the danish north sea,” Ph.D. dissertation, Technical University of Denmark, Kgs. Lyngby, Denmark, 2019.
  * [19] J. B. Jørgensen, “Moving horizon estimation and control,” Ph.D. dissertation, Technical University of Denmark, Kgs. Lyngby, Denmark, 2004.
  * [20] J. B. Jørgensen, G. Frison, N. F. Gade-Nielsen, and B. Damman, “Numerical methods for solution of the extended linear quadratic control problem,” _IFAC Proceedings Volumes_ , vol. 45, no. 17, pp. 187–193, 2012, 4th IFAC Conference on Nonlinear Model Predictive Control.
  * [21] J. B. Jørgensen and S. B. Jørgensen, “Comparison of Prediction-Error-Modelling Criteria,” in _American Control Conference 2007_ , 2007, pp. 140–146.
  * [22] ——, “Continuous-Discrete Time Prediction-Error Identification Relevant for Linear Model Predictive Control,” in _European Control Conference 2007_ , 2007, pp. 4752–4758.
  * [23] ——, “MPC-Relevant Prediction-Error Identification,” in _American Control Conference 2007_ , 2007, pp. 128–133.
  * [24] S. Lenhart and J. Workman, _Optimal Control Applied to Biological Models_ , ser. Chapman & Hall/CRC Mathematical and Computational Biology. Taylor & Francis, 2007.
  * [25] L. Ljungqvist and T. J. Sargent, _Recursive Macroeconomic Theory_. The MIT Press, 2012.
  * [26] S. Martin, “Convergence results for the discrete regularization of linear-quadratic control problems with bang–bang solutions,” _Computational Optimization and Applications_ , vol. 61, no. 3, pp. 731–760, 2015.
  * [27] R. H. Middleton and G. C. Goodwin, _Digital Control and Estimation: A Unified Approach_. Prentice Hall Professional Technical Reference, 1990.
  * [28] C. Moler and C. Van Loan, “Nineteen dubious ways to compute the exponential of a matrix,” _SIAM Rev._ , vol. 20, no. 4, p. 801–836, oct 1978.
  * [29] ——, “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later,” _SIAM Review_ , vol. 45, no. 1, pp. 3–49, 2003\.
  * [30] R. Rockafellar and S. Uryasev, “Conditional value-at-risk: Optimization approach,” _Stochastic optimization: algorithms and applications_ , pp. 411–435, 01 2001.
  * [31] M. Salgado, R. Middleton, and G. C. Goodwin, “Connection between continuous and discrete riccati equations with applications to kalman filtering,” _Control Theory and Applications, IEE Proceedings D_ , vol. 135, pp. 28 – 34, 02 1988.
  * [32] P. M. P. Stanislav Uryasev, _Stochastic Optimization: Algorithms and Applications_. Springer, 2001.
  * [33] S. Uryasev, “Conditional value-at-risk: optimization algorithms and applications,” in _Proceedings of the IEEE/IAFE/INFORMS 2000 Conference on Computational Intelligence for Financial Engineering (CIFEr) (Cat. No.00TH8520)_ , 2000, pp. 49–57.
  * [34] C. Van Loan, “Computing integrals involving the matrix exponential,” _IEEE Transactions on Automatic Control_ , vol. 23, no. 3, pp. 395–404, 1978\.
  * [35] F. L. Walter Alt, Robert Baier and M. Gerdts, “Approximations of linear control problems with bang-bang solutions,” _Optimization_ , vol. 62, no. 1, pp. 9–32, 2013.
  * [36] J. I. Yuz and G. C. Goodwin, _Sampled-Data Models for Linear and Nonlinear Systems (Communications and Control Engineering)_. Springer London, 2014.