# On the Efficient Implementation of an Implicit Discrete-Time Differentiator

J. E. Carvajal-Rubio13, J. D. Sánchez-Torres2, M. Defoort3, M. Djemai3 and A.
G. Loukianov1  1Dept. of Electrical Engineering, CINVESTAV-IPN Guadalajara,
Zapopan, México.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>2Dept. of
Mathematics and Physics, ITESO, Tlaquepaque, México.
Email<EMAIL_ADDRESS>3LAMIH, CNRS UMR 8201, Polytechnic University of
Hauts-de-France, Valenciennes, France.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

New methodologies are designed to reduce the time complexity of an implicit
discrete-time differentiator and the simulation time to implement it. They
rely on Horner’s method and the Shaw-Traub algorithm. The algorithms are
compared for differentiators of order 3, 7, and 10. The Half-Horner and Full-
Horner methods showed the best performance and time complexity.

###### Index Terms:

Time complexity, Implicit Discretization, Differentiator, Sliding-Mode

##  I Introduction

An online differentiator is useful for several applications, such as control
laws based on derivatives of a signal, estimation of unmeasured states and
parameters [1, 2, 3]. The well-known homogeneous differentiator [4] was
proposed in continuous-time and allow to estimate the first $n$ derivatives of
a signal, if its $n$-th derivative has a known Lipschitz constant $L>0$. As
the standard differentiator is implemented in digital systems [5], a discrete-
time version is implemented, for instance the explicit and implicit discrete-
time realizations in [6, 7, 8, 9, 10]. Particularly, the implicit discrete-
time differentiator, proposed in [10], has a remarkable reduction of the
numerical chattering and preserves the main properties of the continuous-time
differentiator [4], i.e., homogeneity property, asymptotic accuracy and
convergence of its errors in finite-time to a vicinity of the origin. The main
drawback of the implicit-discrete time differentiator is that the root of a
polynomial has to be calculated almost each iteration for its implementation.
The polynomial can not be calculated previous to its implementation because
one of its parameters is updated each iteration. To implement this
differentiator, in [10] was proposed the Halley’s method [11], which converges
to the unique positive root of the polynomial with $3$-order for a non-
restrictive set of initial conditions.

Although Halley’s method reduces the time required to implement the implicit
differentiator, Halley’s method needs the evaluation of the function and its
first two derivatives. Then Halley’s algorithm has a quadratic time with
respect to $n$. Several algorithm has been proposed to reduce the number of
basic operations (time complexity [12]) required to evaluate a polynomial and
its derivatives, for instance, Horner’s method [13], Shaw-Traub algorithm
[14], and the De Jong Van algorithm [15]. Other algorithms adapt or
precondition parameters [16], they will not be considered because in this work
the polynomials are updated each iteration. On the other hand, the evaluation
of the implicit differentiator, after calculate the respective roots, presents
a cubic time, but it can be reduced to a quadratic time or linearithmic time
using Horner’s and the discrete Fourier transform [17], respectively.

The first contribution of this paper is to introduce new methodologies
designed to implement the implicit discrete-time differentiator [18] (HIDD),
which allow to reduce the time complexity of the differentiator. The second
one is a numeric comparative between them. This work is organized as follows.
In Section II, implicit differentiator is defined. Additionally, Section III
contains the proposed algorithms for the implementation of implicit
differentiator and its time complexity is calculated. Section IV aims to
compare the time complexity of the algorithms and the simulation time for its
implementation. In Section V, the main results of the paper are summarized and
the future work is presented.

##  II Implicit Discrete-time Differentiator

Let $\tau$ be the sampling time constant and defined as $\tau=t_{k+1}-t_{k}$
for $k=0,1,2,\cdots$. Then HIDD is defined by the following equations:

$\displaystyle\begin{split}&\bm{z}_{k+1}=\bm{\Phi}\left(\tau\right)\bm{z}_{k}+\bm{B}^{\ast}\left(\tau\right)\bm{v}\left(\widetilde{\sigma}_{0,k+1}\right),\\\
&\bm{v}\left(\widetilde{\sigma}_{0,k+1}\right)=\left[\Psi_{0,n}\left(\widetilde{\sigma}_{0,k+1}\right)\;\;\cdots\;\;\Psi_{n,n}\left(\widetilde{\sigma}_{0,k+1}\right)\right]^{T},\\\
&\Psi_{i,n}(\widetilde{\sigma}_{0,k+1})=-\lambda_{n-i}L^{\frac{i+1}{n+1}}\left|\widetilde{\sigma}_{0,k+1}\right|^{\frac{n-i}{n+1}}\xi_{k}.\end{split}$
(1)

For HIDD, $z_{i,k+1}$ corresponds to the estimation of $f^{i}_{0,k}$. On the
other hand, $\bm{\Phi}\left(\tau\right)$ and $\bm{B}^{\ast}\left(\tau\right)$
have the following representation:

$\displaystyle\begin{split}\bm{\Phi}\left(\tau\right)=\left[\begin{array}[]{cccccc}1&\tau&\frac{\tau^{2}}{2!}&\cdots&\frac{\tau^{n-1}}{\left(n-1\right)!}&\frac{\tau^{n}}{n!}\\\
0&1&\tau&\cdots&\frac{\tau^{n-2}}{\left(n-2\right)!}&\frac{\tau^{n-1}}{\left(n-1\right)!}\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&0&\cdots&1&\tau\\\
0&0&0&\cdots&0&1\end{array}\right],\end{split}$ (2)
$\displaystyle\begin{split}\bm{B^{\ast}}\left(\tau\right)=\left[\begin{array}[]{cccccc}\tau&\frac{\tau^{2}}{2!}&\frac{\tau^{3}}{3!}&\cdots&\frac{\tau^{n}}{n!}&\frac{\tau^{n+1}}{\left(n+1\right)!}\\\
0&\tau&\frac{\tau^{2}}{2!}&\cdots&\frac{\tau^{n-1}}{\left(n-1\right)!}&\frac{\tau^{n}}{n!}\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&\cdots&\tau&\frac{\tau^{2}}{2!}\\\
0&0&0&\cdots&0&\tau\end{array}\right].\end{split}$ (3)

Regarding $\widetilde{\sigma}_{0,k+1}$ and $\xi_{k}$, they are calculated
according to the following lemma:

###### Lemma 1

([10]) Let $a_{l}$ and $b_{k}$ be defined as:

$\displaystyle\begin{split}a_{l}&=\frac{\tau^{n-l+1}}{\left(n-l+1\right)!}\lambda_{l}L^{\frac{n-l+1}{n+1}},\;\;l=0,\ldots,n;\\\
b_{k}&=-\sigma_{0,k}-\sum_{l=1}^{n}\frac{\tau^{l}}{l!}z_{l,k}.\end{split}$ (4)

Then $\widetilde{\sigma}_{0,k+1}\in\mathbb{R}$ and $\xi_{k}$ is the unique
pair $\left(\widetilde{\sigma}_{0,k+1},\xi_{k}\right)$ defined conforming to
the following $3$ cases:

  * •

If $b_{k}>a_{0}$, then $\xi_{k}=\left\\{-1\right\\}$ and
$\widetilde{\sigma}_{0,k+1}=-\left(r_{0}\right)^{n+1}\in\mathbb{R}^{-}$ where
$r_{0}$ is the unique positive root of the following polynomial:

$\displaystyle\begin{split}p\left(r\right)=r^{n+1}+a_{n}r^{n}+\cdots+a_{1}r+\left(-b_{k}+a_{0}\right).\end{split}$
(5)

  * •

If $b_{k}\in[-a_{0},a_{0}]$, then $\widetilde{\sigma}_{0,k+1}=0$ and
$\xi_{k}=\left\\{-\frac{b_{k}}{a_{0}}\right\\}$.

  * •

If $b_{k}<-a_{0}$, then $\xi_{k}=\left\\{1\right\\}$ and
$\widetilde{\sigma}_{0,k+1}=r_{0}^{n+1}\in\mathbb{R}^{+}$ where $r_{0}$ is the
unique positive root of the following polynomial:

$\displaystyle\begin{split}p\left(r\right)=r^{n+1}+a_{n}r^{n}+\cdots+a_{1}r+\left(b_{k}+a_{0}\right).\end{split}$
(6)

Note that under the assumption of a constant sampling time $\tau$, $a_{l}$,
and the elements of $\bm{\Phi}\left(\tau\right)$ and
$\bm{B}^{\ast}\left(\tau\right)$ can be calculated previous to the
implementation of HIDD. Opposite to the explicit discrete-time
differentiators, HIDD requires estimating roots of different polynomials each
iteration. They are estimated by using the Halley’s method, which implies an
evaluation of the polynomials (5), (6) and its first two derivatives multiple
times.

###  II-A Halley’s Method

As it was mentioned previously, HIDD requires estimating $r_{0}$. In [10] the
Halley’s method was proposed as a good alternative, which experimentally
needed $3$ iteration to estimate $r_{0}$. It is applied as follows:

$\displaystyle\begin{split}r_{0,j+1}=r_{0,j}-\frac{2\frac{dp(r)}{dr}|_{r=r_{0,j}}p(r_{0,j})}{2\left(\frac{dp(r)}{dr}|_{r=r_{0,j}}\right)^{2}-\frac{d^{2}p(r)}{dr^{2}}|_{r=r_{0,j}}p(r_{0,j})}.\end{split}$
(7)

As $p\left(\left(b_{k}-a_{0}\right)^{\frac{1}{n+1}}\right)>0$ and
$p\left(\left(-b_{k}-a_{0}\right)^{\frac{1}{n+1}}\right)>0$ for the
polynomials (5) and (6), then
$r_{0}\in\left[0,\left(b_{k}-a_{0}\right)^{\frac{1}{n+1}}\right]$ and
$r_{0}\in\left[0,\left(-b_{k}-a_{0}\right)^{\frac{1}{n+1}}\right]$,
respectively. The above matches with the Cauchy’s bound for the roots of a
polynomial [13]. In [10] was demonstrated that the Halley’s method converges
monotonically to $r_{0}$ with an convergence order $3$, [11], for any initial
$r_{0,0}$ belonging to the previous sets of initial conditions. Hence, the
following initial conditions were used:

$\displaystyle\begin{split}r_{0,0}&=\left(\frac{b_{k}-a_{0}}{2}\right)^{1/(n+1)},\;\;\;\;\;\;\;\;\text{for}\;\;b_{k}>a_{0},\\\
r_{0,0}&=\left(\frac{-b_{k}-a_{0}}{2}\right)^{1/(n+1)},\;\;\;\;\;\text{for}\;\;b_{k}<-a_{0}.\end{split}$
(8)

##  III Main Results

As it was mentioned previously, the objective of this work is to reduce the
time complexity of HIDD. First, the variables $\phi_{i}$ and $\bar{b}^{\ast}$
are defined as:

$\displaystyle\begin{split}\phi_{i}&=\frac{\tau^{i-1}}{(i-1)!},\\\
\bar{b}^{\ast}_{i,j}&=\frac{\tau^{j+1-i}}{(j+1-i)!}\lambda_{n-j+1}L^{\frac{j}{n+1}},\end{split}$
(9)

for $i=1,2,\cdots,n+1$ and $j=i,i+1,i+2,\cdots,n+1$. It allows to rewrite
(LABEL:eq:HIDD) as follows:

  * •

If $b_{k}>a_{0}$,

$\displaystyle\begin{split}z_{i,k+1}&=\sum_{j=i}^{n}\phi_{j-i+1}z_{j,k}+\bar{b}^{\ast}_{i+1,j+1}r_{0}^{n-j},\\\
i&=0,1,\cdots,n.\end{split}$ (10)

  * •

If $b_{k}\in[-a_{0},a_{0}]$,

$\displaystyle\begin{split}z_{0,k+1}&=b_{k}+\sum_{j=0}^{n}\frac{\tau^{j}}{j!}z_{j,k},\\\
z_{i,k+1}&=\bar{b}^{\ast}_{i+1,n+1}\left(\frac{b_{k}}{a_{0}}\right)+\sum_{j=i}^{n}\phi_{j-i+1}z_{j,k},\\\
i&=1,2,\cdots,n.\end{split}$ (11)

  * •

If $b_{k}<-a_{0}$,

$\displaystyle\begin{split}z_{i,k+1}&=\sum_{j=i}^{n}\phi_{j-i+1}z_{j,k}-\bar{b}^{\ast}_{i+1,j+1}r_{0}^{n-j},\\\
i&=0,1,\cdots,n.\end{split}$ (12)

###  III-A Direct Evaluation

The number of additions and subtraction, $N_{A}$ and the number of
multiplications and divisions, $N_{M}$, needed to evaluate $z_{i,k+1}$
directly, after obtain $r_{0}$, are calculated as:

$\displaystyle\begin{split}&N_{A1}(n)=(n+1)^{2},\\\
&N_{M1}(n)=\frac{n^{3}}{6}+n^{2}-\frac{1}{6}n-1.\end{split}$ (13)

Therefore, taking into account the $(n+1)$ assignations of $z_{i,k+1}$, the
time complexity is given as:

$\displaystyle\begin{split}T_{1}(n)=\frac{n^{3}}{6}+2n^{2}+\frac{17}{6}n+1,\end{split}$
(14)

which is a cubic time. On the other hand, the Halley’s method is used
recursively each iteration. To evaluate the derivatives and its derivatives,
one could storage the following variables to reduce the number of operations.

$\displaystyle\begin{split}c_{n+1}&=n+1,\\\
c_{i}&=ia_{i},\;\;\;\;\;\;\;\;\;\;\;\;\text{for}\;\;i=1,2,\cdots,n;\\\
d_{n+1}&=n(n+1),\\\ d_{i}&=i(i-1)a_{i},\;\;\text{for}\;\;i=2,3,\cdots,n.\\\
\end{split}$ (15)

Additionally, $\bar{j}$ is defined as the number of iteration used to estimate
$r_{0}$, then the number of additions, subtraction, multiplications and
divisions used to evaluate the polynomials and its derivatives are given as:

$\displaystyle\begin{split}N_{A2}(n)&=\bar{j}\left(3n+1\right),\\\
N_{M2}(n)&=\bar{j}\left(\frac{3}{2}n^{2}+\frac{3}{2}n\right).\end{split}$ (16)

Since Halley’s method is implemented 3 times each iteration for HIDD,
$\bar{j}=3$. Therefore, taking into account the assignation of values, the
evaluation of (7), its initialization and comparatives:

$\displaystyle\begin{split}T_{2}(n)=\frac{9}{2}n^{2}+\frac{27}{2}n+46.\end{split}$
(17)

where one of the operations is a $(n+1)$-th root and a for-loop was
considered. Hence, the complexity of the algorithm is cubic and it is defined
as:

$\displaystyle\begin{split}T(n)=\frac{n^{3}}{6}+\frac{13}{2}n^{2}+\frac{110}{6}n+48.\end{split}$
(18)

where the multiplications and subtraction needed to evaluate $b_{k}$ were
taking into account.

###  III-B Horner Method

Although the variables $\phi_{\cdot}$, $\bar{b}^{\ast}_{\cdot,\cdot}$,
$c_{\cdot}$, $d_{\cdot}$ reduces the number of basic operations, it does not
reduce time complexity of the realization (LABEL:eq:HIDD) with respect to $n$.
Based on the Horner methodology, one could calculate $z_{i,k}$ as follows:

  * •

If $b_{k}>a_{0}$

$\displaystyle\begin{split}&z_{i,k+1}=\sum_{j=i}^{n}\phi_{j-i+1}z_{j,k}+\ldots\\\
&\ldots+(\cdots((\bar{b}^{\ast}_{i+1,i+1})r_{0}+\bar{b}^{\ast}_{i+1,i+2})\cdots)r_{0}+\bar{b}^{\ast}_{i+1,n+1}.\\\
&i=0,1,\cdots,n.\end{split}$ (19)

  * •

If $b_{k}<-a_{0}$

$\displaystyle\begin{split}&z_{i,k+1}=\sum_{j=i}^{n}\phi_{j-i+1}z_{j,k}-\ldots\\\
&\ldots-(\cdots((\bar{b}^{\ast}_{i+1,i+1})r_{0}+\bar{b}^{\ast}_{i+1,i+2})\cdots)r_{0}+\bar{b}^{\ast}_{i+1,n+1}.\\\
&i=0,1,\cdots,n.\end{split}$ (20)

This methodology presents the following number of basic operations:

$\displaystyle\begin{split}N_{A3}(n)&=(n+1)^{2},\\\
N_{M3}(n)&=n(n+1).\end{split}$ (21)

As $\bm{\Phi}(\tau)$ and $\bm{B}(\tau)$ are Toeplitz matrix [19], the time
complexity of evaluate $z_{i,k+1}$ could be reduce to a linearithmic time
($n\log{n}$) using the discrete Fourier transform. This alternative will be
analyzed in a future work. To evaluate the polynomials and its derivatives,
$n$ is considered greater than $1$. Here two methodologies based on Horner’s
method are analyzed, the first one is evaluate the polynomials and derivatives
as follows:

$\displaystyle\begin{split}p\left(r\right)&=(\cdots((r+a_{n})r+a_{n-1})\cdots)r+a_{0}\pm
b_{k},\\\
\frac{dp(r)}{dr}&=(\cdots((c_{n+1}r+c_{n})r+c_{n-1})\cdots)r+c_{1},\\\
\frac{d^{2}p(r)}{dr^{2}}&=(\cdots((d_{n+1}r+d_{n})r+d_{n-1})\cdots)r+d_{2}.\end{split}$
(22)

The methodology (22) use the following number of basic operations:

$\displaystyle\begin{split}N_{A4}(n)&=\bar{j}\left(3n+1\right),\\\
N_{M4}(n)&=\bar{j}\left(3n-1\right).\end{split}$ (23)

Similar to (17), one obtains:

$\displaystyle\begin{split}T_{4}(n)=18n+43.\end{split}$ (24)

However, one could take advantage of the evaluation of $p(r)$ to evaluate
$\frac{dp(r)}{dr}$ and $\frac{d^{2}p(r)}{dr^{2}}$ with the following
methodology for $n\geq 2$:

$\displaystyle\begin{split}F_{i+1}&=rF_{i}+a_{n-i-1},\\\
dF_{i+1}&=rdF_{i}+F_{i+1},\\\ ddF_{i+1}&=rddF_{i}+dF_{i+1},\end{split}$ (25)

for $i=0,\cdots,n-3$, with $F_{0}=r+a_{n}$, $dF_{0}=r+F_{0}$,
$ddF_{0}=r+dF_{0}$, and the value of the evaluation is given as:

$\displaystyle\begin{split}F_{n-1}&=rF_{n-2}+a_{1},\\\ p(r)&=rF_{n-1}+a_{0}\pm
b_{k},\\\ \frac{dp(r)}{dr}&=rdF_{n-2}+F_{n-1},\\\
\frac{d^{2}p(r)}{dr^{2}}&=2ddF_{n-2}.\end{split}$ (26)

It increases the number of assignations and additions but reduces the
multiplications, therefore, one obtains the following time complexity:

$\displaystyle\begin{split}T_{5}(n)=36n+55.\end{split}$ (27)

Using the methodologies (LABEL:eq:Horner_1_1), (LABEL:eq:Horner_1_2) and (22)
a quadratic time is obtained, which is defined as:

$\displaystyle\begin{split}T(n)=2n^{2}+24n+46.\end{split}$ (28)

If (26) is used instead of (22), the quadratic time is given as:

$\displaystyle\begin{split}T(n)=2n^{2}+42n+58.\end{split}$ (29)

###  III-C Shaw–Traub Algorithm

Similar to the methodology (26), an algorithm designed to calculate the
normalized derivatives, $\frac{1}{i!}\frac{d^{i}p(r)}{dr^{i}}$, was proposed
in [14]. Even, the Horner method used to evaluate a polynomial is a special
case of this algorithm. It allows to reduce the number of multiplication but
add divisions, assignations and additions. In this work a modified algorithm
is used, which relies on Shaw–Traub algorithm:

$\displaystyle\begin{split}&t_{1}=r,\;\;t_{i}=t_{i-1}r,\;\;\text{for}\;\;i=2,3,\cdots,n;\\\
&T_{i}^{-1}=a_{n-i}t_{n-i},\;\;\text{for}\;\;i=0,1,\cdots,n-1;\\\
&T_{n}^{-1}=a_{0}\pm b_{k},\\\
&T_{0}^{0}=t_{n}r,\;\;T_{1}^{1}=T_{0}^{0},\;\;T_{2}^{2}=T_{0}^{0};\\\
&T_{i}^{j}=T_{i-1}^{j-1}+T_{i-1}^{j},\;\;\text{for}\;\;j=0,1,2;\;i=j+1,\cdots,n+1;\\\
&p(r)=T_{n+1}^{0},\;\;\frac{dp(r)}{dr}=\frac{T^{1}_{n+1}}{t_{1}},\;\;\frac{d^{2}p(r)}{dr^{2}}=2\frac{T^{2}_{n+1}}{t_{2}}.\end{split}$
(30)

Since $T_{0}^{0}=T_{1}^{1}=T_{2}^{2}$, two assignations could be avoided if
$T_{0}^{0}$ is used instead of $T_{1}^{1}$ and $T_{2}^{2}$. Then the algorithm
(LABEL:eq:Shaw) presents an linear time given as:

$\displaystyle\begin{split}T_{6}(n)=30n+70.\end{split}$ (31)

Applying the algorithms (LABEL:eq:Horner_1_1), (LABEL:eq:Horner_1_2) and
(LABEL:eq:Shaw), HIDD presents a quadratic time, which is given as:

$\displaystyle\begin{split}T(n)=2n^{2}+36n+74.\end{split}$ (32)

It is important to mention that in [20], the best parameters for the family of
algorithms presented in [14] were founded. Then the time complexity (32) could
be reduced tuning its parameters but not its order. As Halley’s method does
not use higher-order derivatives, this work will not consider the algorithm
presented in [15], which improve the algorithm proposed in [14] for the $n+1$
normalized derivatives.

##  IV Simulation

In this section, two comparatives are presented. The first one correspond to a
graphical comparative of the required basic operations of the algorithms. The
second one is a comparative of the time of simulation for the 4 methodologies
proposed in this work. The methodology composed of (7), (10), (11), (12) and a
direct evaluation of the polynomials is referenced as direct evaluation, the
methodology composed of (7), (11) (LABEL:eq:Horner_1_1), (LABEL:eq:Horner_1_2)
and (22) is referenced as half-Horner algorithm, the methodology composed of
(7), (11), (LABEL:eq:Horner_1_1), (LABEL:eq:Horner_1_2) and (26) is referenced
as Full-Horner algorithm and (7), (11), (LABEL:eq:Horner_1_1),
(LABEL:eq:Horner_1_2) and (LABEL:eq:Shaw) is referenced as Shaw-Traub
algorithm.

###  IV-A Simulation I

In this simulation the time complexity of the methodologies are evaluated for
$2\leq n\leq 30$. The results are presented in Figure 1. One can see that the
algorithm with less basic operations is Half-Horner algorithm. Furthermore
Shaw-Traub and Full-Horner algorithm have less basic operations than direct
evaluation for $n>4$. It is important note the difference of Direct evaluation
with respect to the other algorithms for a value of $n>=7$, $240$ or more
basic operations with respect to Half-Horner algorithm.

Figure 1: Comparative between the time complexity of the methodologies.

###  IV-B Simulation II

Since the methodologies have a different proportion of basic operations, this
simulation aims to compare the time require to simulate the algorithms.
Additionally, an simulation without the parameters defined in Equations (9)
and (15) is simulated, which is the same methodology than the direct
evaluation without the parameters $\phi_{i}$, $\bar{b}^{\ast}_{i,j}$, $c_{i}$
and $d_{i}$. Here $n=3$, $n=7$ and $n=10$ are considered, the values of the
signal with noise and the constants defined in (9) and (15) are calculated
previous to the simulation. The sampling time is selected as $\tau=0.001\sec$
and $t=2000,\;10000,\;25000,\;50000\sec$. The results are presented in Tables
I-III. The most efficient methods with respect to the simulation time were the
Half-Horner and Full-Horner algorithms, both present a similar performance for
$n=3$, $n=7$ and $n=10$. For $n=3$, direct evaluation has a better performance
than Shaw-Traub, it contrasts to the results obtained for $n=10$, where Shaw-
Traub algorithm reduces the simulation time compared to direct evaluation. The
above fact matches with its time complexity.

###### Remark 1

Half-Horner and Full-Horner algorithms reduced the simulation time more than
25 times for $n=10$. It can be seen in Table III. It comes from the use of the
variables $\bar{b}^{\ast}_{i,j}$, $\phi_{i}$, $c_{i}$ and $d_{i}$ and a
reduction of the time complexity.

| $2000\sec$ | $10000\sec$ | $25000\sec$ | $50000\sec$ |
---|---|---|---|---|---
Evaluation |  |  |  |  |
without $\phi_{i}$, | $0.5963\sec$ | $2.990\sec$ | $7.470\sec$ | $15.059\sec$ |
$\bar{b}^{\ast}_{i,j}$, $c_{i}$ and $d_{i}$ |  |  |  |  |
Direct |  |  |  |  |
Evaluation. | $0.3578\sec$ | $1.783\sec$ | $4.475\sec$ | $9.027\sec$ |
Half-Horner. | $0.3494\sec$ | $1.753\sec$ | $4.411\sec$ | $8.896\sec$ |
Full-Horner. | $0.3496\sec$ | $1.756\sec$ | $4.407\sec$ | $8.908\sec$ |
Shaw-Traub. | $0.3852\sec$ | $1.922\sec$ | $4.813\sec$ | $9.661\sec$ |
TABLE I: Simulation time of the algorithms for $n=3$ and $\tau=0.001\sec$. | $2000\sec$ | $10000\sec$ | $25000\sec$ | $50000\sec$ |
---|---|---|---|---|---
Evaluation |  |  |  |  |
without $\phi_{i}$, | $6.791\sec$ | $33.808\sec$ | $85.045\sec$ | $169.95\sec$ |
$\bar{b}^{\ast}_{i,j}$, $c_{i}$ and $d_{i}$ |  |  |  |  |
Direct |  |  |  |  |
Evaluation. | $0.486\sec$ | $2.414\sec$ | $6.035\sec$ | $12.51\sec$ |
Half-Horner. | $0.466\sec$ | $2.286\sec$ | $5.729\sec$ | $11.61\sec$ |
Full-Horner. | $0.457\sec$ | $2.293\sec$ | $5.763\sec$ | $11.51\sec$ |
Shaw-Traub. | $0.503\sec$ | $2.46\sec$ | $6.210\sec$ | $12.718\sec$ |
TABLE II: Simulation time of the algorithms for $n=7$ and $\tau=0.001\sec$. | $2000\sec$ | $10000\sec$ | $25000\sec$ | $50000\sec$ |
---|---|---|---|---|---
Evaluation |  |  |  |  |
without $\phi_{i}$, | $14.312\sec$ | $71.6\sec$ | $179.37\sec$ | $358.09\sec$ |
$\bar{b}^{\ast}_{i,j}$, $c_{i}$ and $d_{i}$ |  |  |  |  |
Direct |  |  |  |  |
Evaluation. | $0.831\sec$ | $4.192\sec$ | $10.33\sec$ | $20.527\sec$ |
Half-Horner. | $0.5437\sec$ | $2.75\sec$ | $6.858\sec$ | $13.692\sec$ |
Full-Horner. | $0.5631\sec$ | $2.807\sec$ | $6.997\sec$ | $14.139\sec$ |
Shaw-Traub. | $0.6101\sec$ | $3.097\sec$ | $7.767\sec$ | $15.464\sec$ |
TABLE III: Simulation time of the algorithms for $n=10$ and $\tau=0.001\sec$.

##  V Conclusion

In this work four methodologies were designed to implement the implicit
discrete-time differentiator HIDD, which rely on the Horner’s method and the
Shaw-Traub algorithm. 3 methodologies allow to reduce the time complexity of
the implicit differentiator with respect to $n$, i.e., they present a
quadratic time instead of a cubic time. The methodologies show a noticeable
performance compared to an direct implementation without the parameters
$\phi_{i}$, $\bar{b}_{i,j}$, $c_{i}$ and $d_{i}$. However, the simulations
show that Half-Horner and Full-Horner algorithms reduce the number of basic
operations and its simulation time. Even, they reduced the simulation time of
the differentiator more than 25 times for $n=10$. As future work, the discrete
Fourier transform and the methodology proposed in [15] will be considered to
reduce the time complexity.

## References

  * [1] P. Kaveh and Y. B. Shtessel, “Blood glucose regulation using higher-order sliding mode control,” _International Journal of Robust and Nonlinear Control_ , vol. 18, no. 4‐5, pp. 557–569, 2008.
  * [2] Y. B. Shtessel, I. A. Shkolnikov, and A. Levant, “Smooth second-order sliding modes: Missile guidance application,” _Automatica_ , vol. 43, no. 8, pp. 1470–1476, 2007.
  * [3] M. Iqbal, A. I. Bhatti, S. I. Ayubi, and Q. Khan, “Robust parameter estimation of nonlinear systems using sliding-mode differentiator observer,” _IEEE Transactions on Industrial Electronics_ , vol. 58, no. 2, pp. 680–689, Feb 2011.
  * [4] A. Levant, “Higher-order sliding modes, differentiation and output-feedback control,” _International Journal of Control_ , vol. 76, no. 9-10, pp. 924–941, 2003.
  * [5] V. I. Utkin, J. Guldner, and J. Shi, _Sliding Mode Control in Electro-Mechanical Systems_ , 2nd ed. CRC Press, 2009.
  * [6] M. Livne and A. Levant, “Proper discretization of homogeneous differentiators,” _Automatica_ , vol. 50, no. 8, pp. 2007–2014, 2014.
  * [7] S. Koch and M. Reichhartinger, “Discrete-time equivalent homogeneous differentiators,” in _2018 15th International Workshop on Variable Structure Systems (VSS)_ , July 2018, pp. 354–359.
  * [8] S. Koch, M. Reichhartinger, M. Horn, and L. Fridman, “Discrete-time implementation of homogeneous differentiators,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 2, pp. 757–762, Feb 2020.
  * [9] J.-P. Barbot, A. Levant, M. Livne, and D. Lunz, “Discrete differentiators based on sliding modes,” _Automatica_ , vol. 112, p. 108633, 2020.
  * [10] J. E. Carvajal-Rubio, J. D. Sánchez-Torres, M. Defoort, M. Djemai, and A. G. Loukianov, “Implicit and explicit discrete-time realizations of homogeneous differentiators,” _International Journal of Robust and Nonlinear Control_ , 2021, Special Issue on Homogeneous Sliding-Mode Control and Observation. [Online]. Available: https://doi.org/10.1002/rnc.5505
  * [11] J. M. McNamee and V. Pan, _Numerical Methods for Roots of Polynomials - Part II_ , ser. Studies in Computational Mathematics. Amsterdam London: Elsevier Science, 2013, vol. 16.
  * [12] M. Sipser, _Introduction to the Theory of Computation_. Cengage learning, 2012.
  * [13] J. M. McNamee, _Numerical Methods for Roots of Polynomials - Part I_ , ser. Studies in Computational Mathematics. Amsterdam London: Elsevier Science, 2007, vol. 14.
  * [14] M. Shaw and J. F. Traub, “On the number of multiplications for the evaluation of a polynomial and some of its derivatives,” _Journal of the ACM (JACM)_ , vol. 21, no. 1, pp. 161–167, 1974.
  * [15] L. De Jong and J. Van Leeuwen, “An improved bound on the number of multiplications and divisions necessary to evaluate a polynomial and all its derivatives,” _ACM SIGACT News_ , vol. 7, no. 3, pp. 32–34, 1975.
  * [16] D. E. Knuth, _Art of computer programming, volume 2: Seminumerical algorithms_. Addison-Wesley Professional, 2014.
  * [17] G. H. Golub and C. F. Van Loan, _Matrix computations_ , 4th ed. The Johns Hopkins University Press, Baltimore, 2013\.
  * [18] J. E. Carvajal-Rubio, A. G. Loukianov, J. D. Sánchez-Torres, and M. Defoort, “On the discretization of a class of homogeneous differentiators,” in _2019 16th International Conference on Electrical Engineering, Computing Science and Automatic Control (CCE)_ , September 2019, pp. 1–6.
  * [19] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, _Templates for the solution of algebraic eigenvalue problems: a practical guide_. Society for Industrial and Applied Mathematics (SIAM), 2000.
  * [20] M. Shaw and J. F. Traub, “Analysis of a family of algorithms for the evaluation of a polynomial and some of its derivatives,” _Department of Computer Science Carnegie Mellon University_ , 1975.


# On the improvement of model-predictive controllers.

1st Leander J. Féret  Think Tank
JUMO GmbH & Co. KG
Fulda, Germany
<EMAIL_ADDRESS>2nd Alexander Gepperth  Applied Computer Science
Department
Fulda University of Applied Sciences
Fulda, Germany
0000-0003-2216-7808  3rd Steven Lambeck Fulda University of Applied Sciences
Fulda, Germany
<EMAIL_ADDRESS>Electrical Engineering Department

###### Abstract

This article investigates synthetic model-predictive control (MPC) problems to
demonstrate that an increased precision of the internal prediction model (PM)
automatially entails an improvement of the controller as a whole. In contrast
to reinforcement learning (RL), MPC uses the PM to predict subsequent states
of the controlled system (CS), instead of directly recommending suitable
actions. To assess how the precision of the PM translates into the quality of
the model-predictive controller, we compare a DNN-based PM to the optimal
baseline PM for three well-known control problems of varying complexity. The
baseline PM achieves perfect accuracy by accessing the simulation of the CS
itself. Based on the obtained results, we argue that an improvement of the PM
will always improve the controller as a whole, without considering the impact
of other components such as action selection (which, in this article, relies
on evolutionary optimization).

###### Index Terms:

MPC, model-predictive control, machine learning, DNN, deep-neural-network

##  I Introduction

The use of machine learning (ML) is widespread in different scientific and
industrial use-cases. One of these areas is control science, in which data-
driven ideas have long been prevalent, for example in data-driven system
identification [5], [6].

Even some of the best known use-cases of ML have their roots in control
science, as many will associate ML with self-trained robots learning to walk
[9], [10], [11] or fully automated industry productions [12]. Therefore, ML is
not a new topic in the field of control science.

In general, when systems are difficult to control with a traditional PID-
controller or when the creation of a mathematical model for an MPC [1] or the
training of reinforcement learning (RL) is problematic, the use of data-driven
system-identification is a common option.

With the recent advances in neural networks, many researchers are now working
on using neural networks as PMs for ML-MPCs [2], [3], [7], [8].

This article provides support for the claim that the quality of ML-based
model-predictive control mainly (ML-MPC) depends on the prediction accuracy of
its PM. From this, it is argued, that current and upcoming research on this
topic, can be strictly focused on the core components of an MPC without
showing the improvements of the whole in practical setups.

To provide a strong base for this argument, we employ a basic ML-MPC
architecture as shown in Fig. 1.

Figure 1: A simplified control loop for an MPC and its core components.

This allows us to eliminate any possible architectural reasons for the loss of
control quality, except for the differences resulting from the different PMs.

Consequently, only two points of interest remain, which are further discussed
in section V and VII.

For an overview of the experimental setup, all used core parts of both ML-MPC
and comparison MPC (C-MPC), are discussed in section II, III, IV and V.
Results are presented in VI.

Within VII and VIII it is stated that the prediction accuracy is the main
reason for any given control discrepancy between the ML-MPC and the C-MPC.

###  I-A Related Work

As this article shows, ML-MPC can be split into data-driven system-
identification, non-linear optimization and model-predictive control. Thus,
every contribution revolving around these topics is directly related to the
building of an ML-MPC.

For example, data-driven system-identification is discussed in [13], [14], and
non-linear optimization in [15], and within data-science in [16], [19], [20].

MPC is a stable technology within control engineering, while ML-MPC is
discussed in several domains [18], [21]. It is often analyzed regarding
different prediction models [21], [22], [23]. A prominent discussion of this
aspect is given in [17], which advocates the usage of an LSTM system
customized with dropout layers to better predict noisy data, while preventing
overfitting and using it for model-predictive control. It is shown that
dropout LSTM and co-teaching LSTM give better results on non-Gaussian noisy
data than the standard LSTM network does. While the core topic is data-driven
system identification, the trained PMs are still used in the context of ML-
MPC, and applied to a real chemical system. This is technically unnecessary,
as, as shown in this article, any proven better network prediction is also
directly improving any ML-MPC. The experimental setup could just be used to
show how the prediction improved, without the need for an ML-MPC. This would
reduce information overhead and complexity.

###  I-B Contribution

First of all, this article aims to give a good overview of the basic
architecture of a simple ML-MPC and C-MPC and its components.

More importantly, it suggests that the problem of finding good ML-MPCs can be
split into sub-problems for future research, while giving an experimental
evidence to justify this decision.

This way, the next iteration of work on ML-MPC can focus on the issues ahead,
like the research on better system-identification or the work on non-linear
optimization algorithms. This will reduce the overhead of work and complexity,
and will hopefully result in cleaner and more minimalistic future problem
formulations within this field of research.

##  II Simulations

In classical control theory, non-linear systems are often controlled by
linearizing around fixed points and then applying solutions from linear
control theory. This type of linearization can become quite difficult where CS
are strongly non-linear. Therefore, ML may be able to generate easier and even
better control solutions. To generate training data and targets for the DNN,
the following three simulations of non-linear physical dynamic systems with
varying complexity are used. They will then be used as the CS, as well as the
the PM for C-MPC.

###  II-A Pendulum

The first simulation is an idealized, mathematical model of a pendulum,
defined as a bob of constant mass $\mathit{m[kg]}$ on a rigid, weightless rod
of a constant length $\mathit{l[m]}$, which is fixed to a pivot point with one
degree of freedom as shown in Fig. 2.

Figure 2: Idealized, mathematical pendulum.

Here, the gravitational force is defined as $\mathit{F_{g}=m\cdot g}$, where
$\mathit{g}$ is the acceleration of gravity on earth.

The deflection angle $\varphi[^{\circ}]$ is defined as the angle between the
rod and the rest position of the pendulum, which is the position perpendicular
below the pivot point.

If the pendulum is not at the rest position, a restoring force
$\mathit{F_{tan}}$ will act on it, pointing towards the rest position. The
restoring force increases proportional to the deflection angle $\varphi$. As
the force vector is pointing toward the rest position, (1) has a negative
sign.

$F_{tan}(t)=-m\cdot g\cdot\sin(\varphi(t))$ (1)

Since the pendulum can only move in a circle around the pivot point, the
acceleration is tangential, also called a tangential acceleration
$\mathit{a_{tan}}[^{\circ}\cdot s^{-2}]$, which is defined as the angular
acceleration $\ddot{\varphi}[^{\circ}\cdot s^{-2}]$ multiplied by the length
$\mathit{l}$ of the pendulum.

$a_{tan}(t)=l\cdot\ddot{\varphi}$ (2)

According to Newton’s second law, the equation of motion for the restoring
force can be written as:

$F_{tan}(t)=m\cdot a_{tan}(t)$ (3)

Using the previous equations, the relationship can be written as:

$\ddot{\varphi}+\frac{g}{l}\cdot\sin(\varphi(t))=0$ (4)

According to Newtons First Law, wherein the net force acting on an idealized
point mass is the sum of all involved forces, the next angular velocity
$\mathit{\hat{\dot{\varphi}}}$ can be written as:

$\hat{\dot{\varphi}}=F_{vel}+F_{acc}+F_{in}$ (5)

with the velocity force $\mathit{F_{vel}}$, the angular acceleration force
$\mathit{F_{acc}}$ and an arbitrary input force $\mathit{F_{in}}$.

To calculate the angular velocity force $\mathit{F_{vel}}$, (6) is used, where
the current velocity $\mathit{\dot{\varphi}}$ is adjusted in relation to the
friction $\mathit{b}$, and the pendulum constants, pointing to the resting
position and therefore including a negative sign.

$F_{vel}=-\frac{b}{m\cdot l^{2}}\cdot\dot{\varphi}$ (6)

The angular acceleration force $F_{acc}$ is the result of the gravitational
pull acting on the pendulum towards its rest position as discussed in the
previous equations regarding the tangential restoration force
$\mathit{F_{tan}}$. From this, the angular acceleration
$\mathit{\ddot{\varphi}}$ can be calculated, as shown in (4) and in (7) in the
context of the acceleration force $F_{acc}$.

$F_{acc}=-\frac{g}{l}\cdot\sin(\varphi)$ (7)

The external input force $\mathit{F_{in}}$ is an arbitrary force, given by an
input value $\mathit{u}[N]$ and the physical values of the pendulum shown in
(8).

$F_{in}=\frac{1}{m\cdot l^{2}}\cdot u$ (8)

###  II-B Cartpole

The second model is the mathematical cartpole, which represents a cart of
constant mass $\mathit{m_{c}[kg]}$, that has a pole with a constant mass
$\mathit{m_{p}[kg]}$ connected to it. The pole has a constant length
$\mathit{l_{p}[m]}$ and a theoretical balanced position orthogonal to the top
of the cartpole. The pole angle $\mathit{\theta[^{\circ}]}$ is defined
regarding this balanced position. The cartpole has only one degree of freedom,
where its position $\mathit{x}[m]$ is known. This system is shown in Fig. 3.

Figure 3: Abstracted Cartpole.

Four values represent the state of the system $\mathit{S}$, which are the
cart’s position $\mathit{x}$, the cart’s velocity $\mathit{\dot{x}[m\cdot
s^{-1}]}$, the pole angle $\mathit{\theta[^{\circ}]}$ and the angular velocity
of the pole $\mathit{\dot{\theta}[^{\circ}\cdot s^{-1}]}$.

The simulation of the cartpole is performed in discrete time, and therefore an
arbitrary timestep variable $\tau[s]$ is used.

For each time step within the simulation, the next system state must be
calculated, while using the previous system state and the physical values of
the cartpole. In conclusion, the system state is defined as
$S=[x,\dot{x},\theta,\dot{\theta}]$ and, correspondingly, the next system
state as $\hat{S}=[\hat{x},\hat{\dot{x}},\hat{\theta},\hat{\dot{\theta}}]$.
Its components can be calculated using the following equations:

$\displaystyle\hat{x}$ $\displaystyle=x+\tau\cdot\dot{x}$ (9)
$\displaystyle\hat{\dot{x}}$ $\displaystyle=\dot{x}+\tau\cdot\ddot{x}$
$\displaystyle\hat{\theta}$ $\displaystyle=\theta+\tau\cdot\dot{\theta}$
$\displaystyle\hat{\dot{\theta}}$
$\displaystyle=\dot{\theta}+\tau\cdot\ddot{\theta}$

As the system state $\mathit{S}$ is monitored, only the cart acceleration
$\ddot{x}[m\cdot s^{-2}]$ and the pole angle angular acceleration
$\ddot{\theta}[^{\circ}\cdot s^{-2}]$ must be calculated. Both calculations
are further discussed in [4].

As the cartpole system is controlable, an arbitrary input force $F_{cin}[N]$
can be applied, which is pointing in the same direction as the single degree
of freedom of the cart. With the input force $F_{cin}$, the current system
state $\mathit{S}$ and the physical constants of the cartpole, the angular
acceleration $\ddot{\theta}$ of the pole can be computed as:

$\ddot{\theta}=\frac{g\cdot\sin(\theta)-\cos(\theta)\cdot\frac{F_{cin}+(m_{p}\cdot
l_{p})\cdot\dot{\theta}^{2}\cdot\sin(\theta)}{m_{p}\cdot
m_{c}}}{l_{p}\cdot\frac{4}{3}-\frac{m_{p}\cdot\cos(\theta)^{2}}{m_{p}\cdot
m_{c}}}$ (10)

The cart’s acceleration $\mathit{\ddot{x}}$ can be calculated using (11).

$\ddot{x}=\frac{F_{cin}+(m_{p}\cdot
l_{p})\cdot\dot{\theta}^{2}\cdot\sin(\theta)}{m_{p}\cdot
m_{c}}-\frac{(m_{p}\cdot l_{p})\cdot\ddot{\theta}\cdot\cos(\theta)}{m_{p}\cdot
m_{c}}$ (11)

Given the previous equations, the next state $\hat{S}$ can be calculated,
while dismissing any possible friction.

###  II-C Three-Tank

The three-tank system consists of three connected tanks, as shown in Fig. 4.
$Tank_{1}$ and $Tank_{3}$ are connected to an input valve. The level
$x_{2}[m]$ in $Tank_{2}$ can only change by the flow through the pipes
conneting it to $Tank_{1}$ ($q_{12}$) and $Tank_{3}$ ($q_{23}$). The outgoing
flow is only possible through pipe $q_{3}$ connected to $Tank_{3}$.

Figure 4: Abstracted three-tank system.

For each of the tanks, the change of contained mass $\dot{M_{t}}$ is defined
by any incoming flow of mass $F_{tin}$ subtracted by any outgoing flow of mass
$F_{tout}$, as shown in (12).

$\dot{M_{t}}=F_{tin}-F_{tout}$ (12)

The change of mass $\mathit{\dot{M_{t}}}$ depends on the constant cross-
sectional area of the tank $\mathit{A[m^{2}]}$, the density of the fluid
$\mathit{p[kg\cdot m^{-3}]}$ and the change of the fill level
$\mathit{\dot{x}[m\cdot s^{-1}]}$, as shown in (13).

$\dot{M_{t}}=A\cdot p\cdot\dot{x}(t),$ (13)

where the cross-sectional area of the tank $\mathit{A}$ and the density of the
fluid $\mathit{p}$ are considered constant. In such idealized scenarios, the
mass of the fluid is proportional to its volume.

The change of volume $\dot{V}[m^{3}\cdot s^{-1}]$ depends on the cross-
sectional area of the tank $\mathit{A}$ and the change of the level
$\dot{x}[m\cdot s^{-1}]$, which in return depends on the in- and outgoing
volume-flows $q_{in}(t)[m^{3}\cdot s^{-1}]$ and $q_{out}(t)[m^{3}\cdot
s^{-1}]$, as shown in (14).

$\displaystyle\dot{V}$
$\displaystyle=A\cdot\dot{x}=\sum{q_{in}(t)}-\sum{q_{out}(t)}$ (14)
$\displaystyle\dot{x}$
$\displaystyle=\frac{1}{A}\cdot(\sum{q_{in}(t)}-\sum{q_{out}(t)})$

This can be used to calculate the change of fill level $\dot{x}$ of the
individual tanks. Every tank has different in- and outgoing flows, as shown in
(15).

$\displaystyle\dot{x_{1}}$
$\displaystyle=\frac{1}{A_{1}}\cdot(q_{in_{1}}(t)-q_{out_{12}}(t))$ (15)
$\displaystyle\dot{x_{2}}$
$\displaystyle=\frac{1}{A_{2}}\cdot(q_{in_{12}}(t)-q_{out_{23}}(t))$
$\displaystyle\dot{x_{3}}$
$\displaystyle=\frac{1}{A_{3}}\cdot(q_{in_{3}}(t)+q_{in_{23}}(t)-q_{out_{3}}(t))$

To calculate the in- and outgoing volume flows, the connected pipes must be
simulated. As all of those are very short, their flow dynamics will not be
taken into account. Furthermore, the cross section of the pipes is very small
in comparison to the cross section of the tanks. This results in volume flows
through the connecting pipes between the tanks that only depend on the cross
section and the pipe-end pressures $\mathit{P_{A_{i}}}$ and
$\mathit{P_{A_{j}}}$, which are, in this idealized construction, equal to the
pressure at ground level of the tanks.

For a constant cross section, Bernoulli’s equation (pressure equation) can
describe the volume flow through the pipes with the cross section
$\mathit{a_{ij}}$ of the pipe, an outflow constant $\mathit{\alpha_{ij}}$ with
$0<=\alpha_{ij}>=1$, the fluid density $\mathit{p}$ and the ground pressures
$\mathit{P_{A_{i}}}$ and $\mathit{P_{A_{j}}}$ of the tanks.

$\displaystyle sgn(t)$ $\displaystyle=sign(P_{A_{j}}(t)-P_{A_{i}}(t))$ (16)
$\displaystyle q_{ij}(t)$ $\displaystyle=\alpha_{ij}\cdot a_{ij}\cdot
sgn(t)\cdot\sqrt{\frac{2}{p}\cdot P_{A_{i}}(t)-P_{A_{j}}(t)}$

Furthermore, the presssures at ground level $P_{A}$ depend on the current
level $x$ of the connected tanks. The pressures at ground level $P_{A}$ can be
written as the force acting on the ground of the tank $F_{A}[N]$ regarding the
cross sectional area $A$. The acting ground force $F_{A}$ is the product of
the containing masses $M_{t}$ and the gravitational force $g$. As discussed
before, the containing masses $M_{t}$ can also be expressed by a volume $V$,
given the idealized system architecture. The volume $V$ is calculated by using
the cross sectional area $A$ and the level of the tank $x$. Therefore, it is
possible to simplify the fraction by the cross sectional area $A$.

$\displaystyle P_{A}(t)$ $\displaystyle=\frac{F_{A}(t)}{A}=\frac{M_{t}(t)\cdot
g}{A}=\frac{p\cdot V(t)\cdot g}{A}$ (17) $\displaystyle\frac{p\cdot V(t)\cdot
g}{A}$ $\displaystyle=p\cdot\frac{A\cdot x(t)\cdot g}{A}=p\cdot g\cdot x(t)$

Inserting (17) in (16), the complete flow equation can be written as:

$\displaystyle sgn_{ij}(t)$ $\displaystyle=sign(x_{i}(t)-x_{j}(t))$ (18)
$\displaystyle q_{ij}(t)$ $\displaystyle=\alpha_{ij}\cdot a_{ij}\cdot
sgn_{ij}(t)\cdot\sqrt{2\cdot g\cdot|x_{i}(t)-x_{j}(t)|}$

The volume flow $q_{3}$ of pipe 3 is the outgoing flow of the system. Its
pressure on the disconnected end equals the atmospheric pressure of zero. The
final calculations of every volume flow of each pipe are shown in (19).

$\displaystyle sgn_{ij}(t)$ $\displaystyle=sign(x_{i}(t)-x_{j}(t))$ (19)
$\displaystyle q_{12}(t)$ $\displaystyle=\alpha_{12}\cdot a_{12}\cdot
sgn_{12}(t)\cdot\sqrt{2\cdot g\cdot|x_{1}(t)-x_{2}(t)|}$ $\displaystyle
q_{23}(t)$ $\displaystyle=\alpha_{23}\cdot a_{23}\cdot
sgn_{23}(t)\cdot\sqrt{2\cdot g\cdot|x_{2}(t)-x_{3}(t)|}$ $\displaystyle
q_{3}(t)$ $\displaystyle=\alpha_{3}\cdot a_{3}\cdot\sqrt{2\cdot g\cdot
x_{3}(t)}$

Inserting (19) into (15), every change of each individual level $\dot{x}$ can
be calculated.

##  III Data Preperation

The previously described simulations are used to generate training data. The
data format is a list of samples, each of which represents an unique
measurement session of equal length. Each time step contains either the input,
the output, or a combination of both.

Random input is generated for each sample using the configuration given in
table I under section IV. The simulations produce the corresponding output,
and get reset between samples.

Because DNN models require a combination of input and output data for learning
and predicting, the input and output data are merged, as shown in Figure 5.

Figure 5: Data merge example with five time steps.

Now, one time step includes inputs and their corresponding output. But this
would imply the inclusion of targets in the training data. Therefore, input
and output of the merged data are shifted, so that each time step now contains
the current input and the previous output, while the first input and last
output of the sample are deleted, as shown in Fig. 6.

Figure 6: Data shift example with five time steps.

Since the simulated systems are dynamic, the DNN needs information from past
states to be able to predict the future output. This information is termed
lookback and contains the last $L$ previous values. Therefore, data are
grouped into packages, each containing the current time step and its lookback,
as shown in Fig. 7.

Figure 7: Data lookback example with four time steps and a lookback of three.

To get the data into an understandable format for the DNN, each lookback
package gets collapsed into a single tuple, as shown in Fig. 8.

Figure 8: Data lookback collapse example with two time steps.

Since the input data preparation removes some of the recorded time steps, the
corresponding label data must also be ”cleaned”. Therefore, the first $L$ time
steps of the output data must be removed, equal to the number of values $L$ in
each lookback package. This is shown in Fig. 9.

Figure 9: Data label cleanup example for two time steps.

Finally, the samples of training data and targets each get collapsed into a
list of time steps, as shown in Fig. 10.

Figure 10: Data collapse example for three samples.

##  IV Prediction Models

Each MPC requires its own prediction model for the CS under consideration.

The used ML-MPC is using a single DNN architecture for each of the three
problems as its PM, with the configuration shown in table I.

TABLE I: DNN-Model Descriptions General |  Activation: Relu
---|---
Loss: mse
Framework: Keras v. 2.8.0
Random Seed: 0
Pendulum |  Model Summary Layer (type) Output Shape in (InputLayer) [(None, 8)] dense0 (Dense) (None, 16) dense1 (Dense) (None, 12) dense2 (Dense) (None, 8) dense3 (Dense) (None, 4) output (Dense) (None, 1) |  | Configuration
---
Samples |  1056
Steps |  34
Lookback |  3
Test Split |  0.2
Epochs |  64
Normalize |  No
In-/Output Ranges
Force |  -10 - 10
Angle |  -8 - 8
Cartpole |  Model Summary Layer (type) Output Shape in (InputLayer) [(None, 18)] dense0 (Dense) (None, 64) dense1 (Dense) (None, 64) dense2 (Dense) (None, 64) dense3 (Dense) (None, 64) dense4 (Dense) (None, 64) dense5 (Dense) (None, 64) output (Dense) (None, 2) |  | Configuration
---
Samples |  7012
Steps |  12
Lookback |  5
Test Split |  0.2
Epochs |  150
Normalize |  No
In-/Output Ranges
Force |  -20 - 20
Position |  -8 - 8
Angle |  -11 - 11
Tanks |  | Model Summary
---
Layer (type) |  Output Shape
in (InputLayer) |  [(None, 20)]
dense0 (Dense) |  (None, 24)
dense1 (Dense) |  (None, 24)
dense2 (Dense) |  (None, 16)
dense3 (Dense) |  (None, 12)
output (Dense) |  (None, 3)
| Configuration
---
Samples |  806
Steps |  24
Lookback |  3
Test Split |  0.2
Epochs |  80
Normalize |  Yes
In-/Output Ranges
Inflow1-2 |  0 - 0.0001
Level1-3 |  0 - 1.0

The C-MPC uses the underlying simulation itself as its PM. C-MPC therefore
have access to predictions that are guaranteed to produce the best possible
result that the ML-MPC could ever reach.

##  V Model Predictive Control

Any MPC can be formulated as an optimization of the best next input (or:
action, the language of RL) using a model of the CS for future predictions,
based on an arbitrary reference value, that defines the state, the CS should
reach, also often called setpoint in control science.

Since no linearization of the discussed non-linear systems is used, the MPC
includes a non-linear optimization algorithm, in particular, a custom-built
genetic algorithm.

Using this approach, the ML-MPC and the C-MPC are not guaranteed to find the
globally optimal action. However, since this deficiency is shared by both
architectures, the comparison remains valid.

The used control loop is straightforward. The CS receives an input (action)
and returns an output. Then, the MPC searches the best subsequent input based
on the set of recorded outputs and the given reference value. The found input
is now again applied to the CS, thus terminating one control iteration.

Within each control iteration, a non-linear optimization problem must be
solved. Depending on the chosen algorithm, this takes time, energy and often
does not guarantee the global minimum. Within the optimization, the PM has to
make multiple future predictions, which means, a slow prediction time of the
DNN or the simulation is one core issue within the non-linear optimization
algorithm. Furthermore, each prediction within the optimization has to be done
from the state the CS is currently in.

As the C-MPC uses the same simulation used for the CS as its PM, it can just
mirror the simulations states one to one, which allows for 100% prediction
accuracy.

The ML-MPC, on the other hand, uses a DNN as its PM, where the inner states
cannot be set to mirror the CS. Therefore, the current state of the CS is
encoded within the input data. This is done by creating the first input state
for each prediction with the past outputs of the CS and discarding any
possible outputs from the DNN. This also means that the DNN can only start to
predict from the correct state when enough output data of the CS is collected
to create the first lookback-package for the input.

The process of encoding the states of the CS into the input data of the DNN is
further called the state correction, and a DNN with corrected state is called
corrected DNN (C-DNN). Shown is a comparison between a prediction session of a
DNN with and without state correction in Fig. 11.

Figure 11: Prediction session of the pendulum system with and without state
correction.

It is evident that the state correction process is effectivly removing the
summing of prediction errors over time. Furthermore, we observe that the state
of the CS can be encoded into the input data of the DNN.

With state correction, it is possible to start each future prediction within
the optimization from the correct CS state. But this can only be done between
each control iteration, hence the predictions within the optimization include
any summed prediction errors.

However, as the state correction is effective, the DNN can be optimized to
only predict the needed length for the optimization; the prediction horizon.
Furthermore, any quality loss of the control cannot be descriped to possible
incorrect state correction.

##  VI Comparison

The final control comparison is run for each of the simulations, with its
result presented in the following subsections.

###  VI-A Comparison-Result: Pendulum

The control session on the pendulum system is shown in Fig. 12.

Figure 12: ML-MPC control versus C-MPC control with three reference jumps on
the pendulum simulation.

###  VI-B Comparison-Result: Cartpole

The control session on the cartpole system is shown in Fig. 13.

Figure 13: ML-MPC control versus C-MPC control with three reference jumps on
the cartpole simulation.

###  VI-C Comparison-Result: Three-Tank

The control session on the three-tank system is shown in Fig. 14.

Figure 14: ML-MPC control versus C-MPC control with three reference jumps on
the three-tank simulation.

##  VII Discussion

The presented results demonstrate that the control works with varying degrees
of accuracy and quality, while the C-MPC has the better result.

Since the MPCs only differ in their PMs, it is reasonable to assume that any
discrepancy between their results comes from the prediction errors within the
optimization.

This is shown by plotting the used predictions within the optimization and
comparing them to the ideal ones, as done in Fig. 15.

Figure 15: Cartpole control run and zoomed in area of interest with PM
predictions.

Since the state correction was shown to be effective, the prediction errors
within the control simulations are the result of the inherent prediction
errors of the trained DNNs.

##  VIII Conclusion

The experiments show that DNNs can be used for ML-MPCs, and that any
diminished control quality regarding the C-MPC is the result of the prediction
error only.

Knowing this, the ML-MPC should be improved within its separated components,
where a better prediction accuracy of the ML-algorithmn directly results in a
better control quality. There is no need to develop and test on ML-MPCs, as
any compatible MPC-architecture can be improved without an ML-MPC in mind.

Hence, separating the issues into categories of machine learning (prediction
accuracy and prediction time), control (MPC-architectures) and non-linear
optimization (prediction optimization algorithm like genetic algorithm or
quadratic programming) is possible and recommended, while at the same time
guaranteeing an improvement of any ML-MPC.

## Acknowledgment

We thank Dr. Johann Letnev and Dr. Ivan Jursic for useful hints about the
manuscript.

## References

  * [1] R. Dittmar, and B.-M. Pfeiffer, “Modellbasierte Prädiktive Regelung”, 2004, ISBN: 9783486594911.
  * [2] Z. Wu, and P. D. Christofides, “Economic machine learning-Based Predictive Control of Nonlinear Systems,” in Mathematics, 2019.
  * [3] J. Luo, “Machine Learning Modeling and Predictive Control of Nonlin- ear Processes Using Noisy Data”, in Chinese Journal of Aeronautics 35.4, Apr. 2022, pp. 14– 44\.
  * [4] R. V. Florian, “Correct equations for the dynamics of the cart-pole system,” in Center for Cognitive and Neural Studies (Coneural), July 2005.
  * [5] I. R. Manchester, M. Revay, and Ruigang Wang, “Contraction-Based Methods for Stable Identification and Robust Machine Learning: a Tutorial”, in 60th IEEE conference on Decision and Control, 2021.
  * [6] E. Vladu, “Using Genetic Algorithms in System Identification”, in Tech. rep. Departement of Electrical Engineering and Information Technology, University of Oradea, 2016.
  * [7] Z. Wu, “Machine learning-based predictive control of nonlinear processes. Part I: Theory”, AIChE Journal 65, 2019.
  * [8] Z. Wu, “machine learning-based predictive control of nonlinear processes. Part II: Computational implementation”, AIChE Journal 65, 2019.
  * [9] J. Morimoto, and K. Doya, “Reinforcement learning of dynamic motor sequence: Learning to stand up”, Intl. Conference on Intelligent Robots and Systems, Oct. 1998.
  * [10] A. L. Kun, and W. T. Miller, “ADAPTIVE DYNAMIC BALANCE OF A BIPED ROBOT USING NEURAL NETWORKS”, Proceedings of IEEE International Conference on Robotics and Automation, Apr. 1996.
  * [11] S. Ha, P Xu, Z. Tan, S. Levine, and J. Tan, “Learning to Walk in the Real World with Minimal Human Effort”, 4th Conference on Robot Learning, 2020.
  * [12] A. S. Rajawat, R. Rawat, K. Barhanpurkar, R. N. Shaw, and A. Gosh, “Robotic process automation with increasing productivity and improving product quality using artificial intelligence and machine learning”, Artificial Intelligence for Future Generation Robotics, 2021, pp. 1-13.
  * [13] M. Forgione, and D. Piga, “Continuous-time system identification with neural networks: MOdel structures and fitting criteria”, European Journal of Control, 2021.
  * [14] S. L. Brunton, and N. Kutz, “Data-Driven Science and Engineering”, Cambridge University Press, 2019, ISBN: 9781108422093.
  * [15] A. Ruszczyński, “Nonlinear Optimization”, Princeton University Press, 2006, ISBN: 9780691119151.
  * [16] G. Joya, M. A. Atencia, and F. Sandoval, “Hopfield neural networks for optimization: study of the different dynamics”, Neurocomputing 43, 2002.
  * [17] J. Luo, “Machine Learning Modeling and Predictive Control of Nonlinear Processes Using Noisy Data”, UCLA Electronic Theses and Dissertations, 2021.
  * [18] A. Draeger, S. Engell, and H. Ranke, “Model predictive control using neural networks”, IEEE Control Systems Magazine, 1995.
  * [19] Y. Pan, and J. Wang, “Model Predictive Control of Unknown Nonlinear Dynamical Systems Based on Recurrent Neural Networks”, IEEE Transactions on Industrial Electronics, 2012.
  * [20] B. M. Åkesson, H. T. Toivonen, “A neural network model predictive controller”, Journal of Process Control, 2006.
  * [21] P. Kittisupakorn, P. Thitiyasook, M.A. Hussain, and W. Daosud, “Neural network based model predictive control for a steel pickling process”, Journal of Process Control, 2008.
  * [22] Z. Wu, D. Rincon, and P. D. Christofides, “Process structure-based recurrent neural network modeling for model predictive control of nonlinear processes”, Journal of Process Control, 2020.
  * [23] S. Piché, J. Keeler, G. Martin, G. Boe, D. Johnson, and M. Gerules, “Neural Network Based Model Predictive Control”, Advances in Neural Information Processing Systems 12 (NIPS), 1999.