# Torus Queen Independence

Kada K Williams

(April 28, 2024)

###### Abstract

Define a queen on $\mathbb{Z}_{n}^{d}$ with admissible moves parallel to
$\mathbf{x}\in\\{-1,0,1\\}^{d}$ at arbitrary length. How many queens can be
placed on $\mathbb{Z}_{n}^{d}$ without any two in conflict? We improve a 1918
result by Pólya on the torus board $\mathbf{Z}_{n}^{2}$ while deducing an
exact answer whenever $n\leq 15$. Moreover, we give an example of
$(n-c)n^{d-2}$ independent queens on $\mathbf{Z}_{n}^{d}$, where $c(2)=4$.

##  1 Introduction

According to Hedetniemi and Hedetniemi [1], significant attention was devoted
to a wide variety of problems that involve chess pieces upon an array of size
$n\times n$. Each figurine, such as a bishop, a rook, or a queen, has a
prescribed set of moves. If a figure can move to a specific field and capture
it, a figurine there is in check. How many of a piece can be placed on a board
such that all the fields are in check? For example, the minimal number of
kings is $\left\lceil\frac{n+2}{3}\right\rceil^{2}$, via grids of step size
$3$.

Let us now examine how many pieces can be placed of pairwise distinct colours.
If no piece is to hold another in check, what is the maximal number of pieces?
This exceeds the minimum in the previous domination scenario, such as that
considered in Barr and Rao [2], because if a field is out of check, a piece
can be placed there, provided the moves are mutual. For example, in the case
of independent rooks, the maximum is obviously $n$. On the $n\times n$ board,
we still have sufficient wiggle room to place $n$ queens [1]. However, Pólya
[3] showed that if $n$ and $6$ are not coprime, the torus-like
$\mathbb{Z}_{n}\times\mathbb{Z}_{n}$ cannot field $n$ queens.

Our goal in this paper is to carefully analyze how many queens can occupy
$\mathbb{Z}_{n}^{2}$ without conflict, before leaping into $d$ dimensions.
First, we provide constructions of at least $n-4$ independent queens on
$\mathbb{Z}_{n}^{2}$, broadening work by Burger, Cockayne, and Mynhardt [4].
Next, we establish that $n-1$ queens are optimal if and only if $n\equiv\pm
2\pmod{12}$. After that, we tackle the problem of [3] on $\mathbb{Z}_{n}^{3}$,
obtaining a modulo $25$ condition for the impossibility of $n^{2}$ independent
queens. Finally, we shall place $n^{d-1}-c(d)n^{d-2}$ queens on
$\mathbb{Z}_{n}^{d}$ without conflict.

##  2 Constructions

######  Theorem 1.

Let $n$ be a positive integer that is coprime to $2\cdot 3$. Then the maximal
number of queens that can be placed on $\mathbb{Z}_{n}^{2}$ without conflict
is $n$.

###### Proof.

We consider the fields $(t,2t)$, where $0\leq t\leq n-1$, one for each column.
These are in pairwise distinct rows, as we assumed that $n|2t$ implies $n|t$.

If we move from field $(x,y)$ diagonally by $(d,-d)$, the value of $x+y$
remains constant. If we move by $(d,d)$, the value of $x-y$ is constant.
Hence, a diagonal is described by the sum or the difference of its
coordinates.

In our example of $(t,2t)$, the sum of coordinates equals $3t$, while the
difference of coordinates equals $t$. Since $n$ is coprime to $3$, these
values are also pairwise distinct. It follows that $n$ queens placed on these
fields are not in conflict. ∎

######  Theorem 2.

2mod4no3 Let $n$ be a positive integer such that $n\equiv\pm 2\pmod{12}$. Then
it is possible to place $n-1$ queens on $\mathbb{Z}_{n}^{2}$ without conflict.

###### Proof.

In this case, $n=2m$ is even, where $m$ is odd and $3$ does not divide it.
Hence, it is auspicious to consider $(t,3t)$, since the values $3t$ are all
distinct.

Inspecting the diagonals, the sums are $4t$ and the differences $2t$. These
attain the same value mod $n$ exactly if the values of $t$ differ by
$\frac{n}{2}=m$. Only for the values $0\leq t\leq m-1$ are the queens on
$(t,3t)$ not in conflict, see Figure 1.

Thus, if we place queens on $(t-1,3t)$ for $m<t\leq 2m-1$, distinct rows and
columns, then the sum values are $4t-1$, the difference values $2t-1$, which
are odd values, all distinct. Further, since odd and even values are distinct,
we have $m+(m-1)=n-1$ queens placed out of conflict. ∎

Figure 1: $9$ queens on $\mathbb{Z}_{10}^{2}$ according to 2mod4no3

It is worth pondering what happens if we work the fields $(t,2t)$. In Figure
2, we see that if $n$ is a multiple of $3$, only when $0\leq
t\leq\frac{n}{3}-1$ are these not in conflict.

$X$ Figure 2: An optimal arrangement of queens on $\mathbb{Z}_{9}^{2}$

######  Lemma 3.

Let $n$ be an odd multiple of $3$. Then it is possible to place $n-4$ queens
on $\mathbb{Z}_{n}^{2}$ without conflict.

###### Proof.

Let $n=3m$, where $m$ is odd. The fields $(t,2t)$ with $0\leq t\leq m-1$ have
sum values $3t$, all possible values in $\mathbb{Z}_{n}$ that are $0\pmod{3}$.
To these, we add the fields $(t,2t+2)$ with $m\leq t\leq 2m-1$, whose sum
values are all possible values $2\pmod{3}$. These are in distinct rows,
because $m$ and $2$ are coprime. Finally, we include $(t,2t+4)$ with $2m\leq
t\leq 3m-5$ for a total of $n-4$ fields. This way, the difference values are
all distinct, as are the rows, and the sums vary modulo $3$. ∎

Figure 3: An optimal arrangement of queens on $\mathbb{Z}_{15}^{2}$

Figures 2 and 3 depict how small values of $n$ allow as many as $n-2$ queens
on $\mathbb{Z}_{n}^{2}$ without conflict. Later, we shall show these to be
optimal.

######  Lemma 4.

Let $n$ be an even positive integer that is coprime to $3$. Then we can place
$n-4$ queens on $\mathbb{Z}_{n}^{2}$ without conflict.

###### Proof.

Along the diagonals of $(t,2t)$, the sum values are $3t$, increasing by $3$ in
every step. Therefore, we place queens on $(t,2t)$ for $0\leq
t\leq\frac{n}{2}-1$ and $(t,2t+3)$ for $\frac{n}{2}\leq t\leq n-4$, they are
not in conflict. The row values $2t$ and $2t+3$ differ modulo $2$. The
difference values $t$ and $t+3$ do not overlap on the ranges we specified.
Moreover, the sum values $3t$ and $3(t+1)$ are pairwise distinct, because $n$
and $3$ are coprime. If $n-4$ queens occupy these fields, there is no
conflict. ∎

Figure 4: $n-2$ queens on $\mathbb{Z}_{n}^{2}$, for $n=8$ and $n=12$

######  Lemma 5.

Let $n$ be a positive integer multiple of $2\cdot 3$. Then it is possible to
place $n-3$ queens on $\mathbb{Z}_{n}^{2}$ without conflict.

###### Proof.

Writing $n=6m$, we place queens on the following fields:

$\displaystyle\\{(t,2t)|0\leq t\leq 2m-1\\}\cup\\{(t,2t+2)|2m\leq t\leq
3m-2\\}\cup$ $\displaystyle\cup$ $\displaystyle\\{(t,2t+1)|3m\leq t\leq
5m-3\\}\cup\\{(t,2t-1)|5m\leq t\leq 6m-1\\}.$

These $n-3$ queens are in distinct rows, with $\frac{n}{2}-1$ even rows and
$\frac{n}{2}-2$ odd rows. Due to the obvious omissions, the difference values
equal $t$ ($0\leq t\leq 2m-1$), $t+2$ ($2m\leq t\leq 3m-2$), $t+1$ ($3m\leq
t\leq 5m-3$), $t-1$ ($5m\leq t\leq 6m-1$), all distinct, leaving out $2m$,
$2m+1$, and $6m-1$. The sum values are also all distinct, covering all
$0\pmod{3}$ values, all $2\pmod{3}$ values other than $-1$, and all
$1\pmod{3}$ values with two exceptions, $-2$ and $-5$. ∎

######  Corollary 6.

For all $n\geq 16$, there can be $n-4$ queens on $\mathbb{Z}_{n}^{2}$ without
conflict.

This construction is not sharp, see Figures 1, 2, 3, and 4.

######  Conjecture 7.

For all $n\geq 3$, there can be $n-2$ queens on $\mathbb{Z}_{n}^{2}$ without
conflict.

##  3 Impossibility proofs

In 1918, Pólya [3] proved that if $n$ is not coprime to $2\cdot 3$, then it is
impossible to place $n$ queens on $\mathbb{Z}_{n}^{2}$. The reasoning involved
can be sharpened as follows.

######  Lemma 8.

Let $n$ be an odd positive integer, and suppose that $n-1$ queens are placed
on $\mathbb{Z}_{n}^{2}$ without conflict. Then one can place one more queen
without conflict.

###### Proof.

Let our queens occupy the fields $(x_{1},y_{1}),\dots,(x_{n-1},y_{n-1})$. It
is known that all the $x_{i}$, $y_{i}$, $x_{i}\pm y_{i}$ values are distinct.
Our objective is to find which value in $\mathbb{Z}_{n}$ is not attained.

By translating the rows and columns, we may suppose that $x_{i}$ and $y_{i}$
take on every value but $0$. Let us then inquire as for whether a queen can be
placed on $(0,0)$ without conflict.

Denote by $a_{+}$ and $a_{-}$ the missing value of $x_{i}+y_{i}$ and
$x_{i}-y_{i}$, respectively. If we can determine these to equal $0$, then the
diagonals of $(0,0)$ are not in check.

It is known that $\sum_{i=1}^{n-1}x_{i}$ and $\sum_{i=1}^{n-1}y_{i}$ equal
$S_{1}=\sum_{i=1}^{n-1}i=\frac{n(n-1)}{2}$. If $n$ is odd, this is divisible
by $n$. Therefore,

$\sum_{i=1}^{n-1}(x_{i}\pm y_{i})\equiv S_{1}-a_{\pm}\pmod{n}$

yields $a_{\pm}\equiv 0\pmod{n}$, and a queen on $(0,0)$ is established
without conflict. ∎

Since placing $n-1$ queens out of check determines the missing column, row,
difference diagonal, and sum diagonal, the congruences in [3] can be improved.

######  Theorem 9.

Let $n$ be a multiple of $3$ or a multiple of $4$. Then it is impossible to
place $n-1$ queens on $\mathbb{Z}_{n}^{2}$ without conflict.

###### Proof.

Let our queens occupy the fields $(x_{1},y_{1}),\dots,(x_{n-1},y_{n-1})$. It
is known that all the $x_{i}$, $y_{i}$, $x_{i}\pm y_{i}$ values are distinct.

As before, we suppose that $x_{i}$ and $y_{i}$ do not attain the value $0$,
and the missing value for $x_{i}\pm y_{i}$ is
$a_{\pm}\equiv\sum_{i=1}^{n-1}i\pmod{n}$. If $n$ is odd, this is congruent to
$0$, while if $n$ is even, this is congruent to $\frac{n}{2}$ modulo $n$.

The sums of $x_{i}^{2}$ and of $y_{i}^{2}$ are
$S_{2}=\sum_{i=1}^{n-1}i^{2}=\frac{(n-1)n(2n-1)}{6}$. Diagonally,

$\sum_{i=1}^{n-1}(x_{i}\pm y_{i})^{2}\equiv S_{2}-a_{\pm}^{2}\pmod{n}.$

We shall make use of the parallelogram identity

$\sum_{i=1}^{n-1}(x_{i}+y_{i})^{2}+\sum_{i=1}^{n-1}(x_{i}-y_{i})^{2}=2\sum_{i=1}^{n-1}x_{i}^{2}+2\sum_{i=1}^{n-1}y_{i}^{2}.$

Plugging in, we obtain $2\left(S_{2}-a_{\pm}^{2}\right)\equiv 4S_{2}\pmod{n}$,
so $2S_{2}\equiv-2a_{\pm}^{2}\pmod{n}$. If $n$ is an odd multiple of $3$, we
have a contradiction, given $a_{\pm}\equiv 0\pmod{n}$ and

$2S_{2}=\frac{(n-1)n(2n-1)}{3}=\frac{n}{3}\cdot(n-1)(2n-1)\equiv\frac{n}{3}\pmod{n}.$

If $n$ is even, a stronger congruence holds. Observe that a shift of $X$ by
$+n$ increments the value $X^{2}$ by $(X+n)^{2}-X^{2}=2nX+n^{2}$, which is
divisible by $2n$. Hence, the sum of squares congruence actually holds modulo
$2n$, yielding $S_{2}\equiv-a_{\pm}^{2}\pmod{n}$.

If $n$ is divisible by $4$, then $a_{\pm}^{2}=\frac{n^{2}}{4}$ vanishes modulo
$n$. However, $S_{2}$ is not divisible by $n$, because $\frac{(n-1)(2n-1)}{6}$
is not an integer. This is a contradiction.

Since either case is contradictory, a placement of $n-1$ queens is impossible.
∎

##  4 Higher dimensions

######  Definition 10.

A queen on $\mathbb{Z}_{n}^{d}$ can move any number of times by
$\mathbf{x}\in\\{-1,0,1\\}^{d}$.

In the problem of conflict-free queens on $\mathbb{Z}_{n}^{d}$, a restriction
to fields where the last coordinate is fixed reveals a problem of conflict-
free queens on $\mathbb{Z}_{n}^{d-1}$. Therefore, the number of queens on
$\mathbb{Z}_{n}^{d}$ is at most $n$ times the number of queens on
$\mathbb{Z}_{n}^{d-1}$. This provides an upper bound of $n^{d-1}$, achievable
for infinitely many values of $n$.

######  Theorem 11.

Let $n$ be a positive integer that has no prime divisor less than $2^{d}$.
Then $n^{d-1}$ many queens can be placed on $\mathbb{Z}_{n}^{d}$ without
conflict, as many as rooks.

###### Proof.

We place queens on the fields

$\\{(t_{1},-2t_{1}+t_{2},-2t_{2}+t_{3},\dots,-2t_{d-1}):t_{1},\dots,t_{d-1}\in\mathbb{Z}_{n}\\}.$

If two queens are in conflict, then their displacement has the form

$(t_{1},-2t_{1}+t_{2},-2t_{2}+t_{3},\dots,-2t_{d-1})=(\epsilon_{1}t_{d},\epsilon_{2}t_{d},\epsilon_{3}t_{d},\dots,\epsilon_{d}t_{d}),$

where $(\epsilon_{1},\dots,\epsilon_{d})\in\\{-1,0,1\\}^{d}$ is nonzero and
$t_{d}\neq 0$. Solving, we see

$t_{1}=\epsilon_{1}t_{d},\quad t_{2}=(\epsilon_{2}+2\epsilon_{1})t_{d},\quad
t_{3}=(\epsilon_{3}+2\epsilon_{2}+4\epsilon_{1})t_{d},\quad\dots$

culminating in

$0=(\epsilon_{d}+2\epsilon_{d-1}+\dots+2^{d-1}\epsilon_{1})t_{d}.$

The bracket contains a nonzero integer of at most $d$ digits in binary, of
size less than $2^{d}$. Provided that $n$ is coprime to it, $t_{d}=0$ in
$\mathbb{Z}_{n}$, and so the conflict is absurd. ∎

Is it impossible to place $n^{d-1}$ queens if $n$ has a prime divisor smaller
than $2^{d}$? Our congruence method is limited, as we could only obtain a
partial answer to this.

######  Theorem 12.

Let $n$ be a positive integer that is coprime to $2\cdot 3$, a multiple of
$5$, and not a multiple of $5^{2}$. Then it is impossible to place $n^{2}$
many queens on $\mathbb{Z}_{n}^{3}$ without conflict.

###### Proof.

Suppose that queens can be placed on

$(x_{1},y_{1},z_{1}),\dots,(x_{n^{2}},y_{n^{2}},z_{n^{2}})\in\mathbb{Z}_{n}^{3}$

without conflict. A queen can move parallel to $(1,0,0)$, $(1,\pm 1,0)$,
$(1,\pm 1,\pm 1)$, and permutations. Moving $(x,y,z)$ along these, the value
of $(y,z)$, $(x\pm y,z)$, and $(x\pm y,x\pm z)$ shall stay constant,
respectively. Hence, across the queens on $(x_{i},y_{i},z_{i})$, the
respective pairs take on distinct values in $\mathbb{Z}_{n}^{2}$, and so every
value.

Let $p$ and $q$ be arbitrary exponents, and consider the sum

$\sum_{i=1}^{n^{2}}x_{i}^{p}y_{i}^{q}.$

Since $(x_{i},y_{i})$ takes on every possible value in $\mathbb{Z}_{n}^{2}$
once, the sum equals

$\sum_{(x,y)\in\mathbb{Z}_{n}^{2}}x^{p}y^{q}=\left(\sum_{x\in\mathbb{Z}_{n}}x^{p}\right)\left(\sum_{y\in\mathbb{Z}_{n}}y^{q}\right).$

The value of such power sums can be computed recursively via the identity

$\sum_{x=0}^{n-1}\frac{x(x-1)\dots(x-p+1)}{p!}=\frac{n(n-1)\dots(n-p)}{(p+1)!}.$

Indeed, on the left-hand side, the polynomial expression in $x$ is the sum of
$x^{p}$ and lower order terms. For our purpose, it suffices to extract how
$S_{p}=\sum_{x\in\mathbb{Z}_{n}}x^{p}$ is a multiple of $n$ for $p=0,1,2,3$,
but for $p=4$, it happens to be $\frac{n}{5}$ modulo $n$, which is not a
multiple of $5$.

Consider the first such sum that does not vanish modulo $n$, given by

$\sum_{i=1}^{n^{2}}x_{i}^{4}y_{i}^{4}=S_{4}\cdot
S_{4}\equiv\left(\frac{n}{5}\right)^{2}\pmod{n}.$

Dropping the indices, we implement a shorthand for the sum, written as $\sum
x^{4}y^{4}$. For the same reason, we have

$\sum(x\pm y)^{4}(x\pm z)^{4}\equiv S_{4}^{2}\pmod{n},$

and the same is true if we use $y^{4}$ in place of $(x\pm y)^{4}$. Thus,

$\sum\left[(x+y)^{4}+(x-y)^{4}\right]\cdot\left[(x+z)^{4}+(x-z)^{4}\right]\equiv
4\cdot S_{4}^{2}\pmod{n},$
$\sum\left[(x+y)^{4}+(x-y)^{4}-2y^{4}\right]\cdot\left[(x+z)^{4}+(x-z)^{4}-2z^{4}\right]\equiv
0\pmod{n}.$

Expanding $(x+y)^{4}+(x-y)^{4}=2(x^{4}+6x^{2}y^{2}+y^{4})$, the latter
equation shows $\sum(x^{4}+6x^{2}y^{2})(x^{4}+6x^{2}z^{2})$ to vanish modulo
$n$. Of course, we know this to be true of $\sum x^{6}z^{2}\equiv
S_{6}S_{2}\equiv 0\pmod{n}$, as well as $\sum x^{8}$ and $\sum x^{6}y^{2}$.
Therefore, we learn that $\sum x^{4}y^{2}z^{2}$ vanishes modulo $n$. Returning
to the former equation,

$\sum(x^{4}+6x^{2}y^{2}+y^{4})(x^{4}+6x^{2}z^{2}+z^{4})\equiv
S_{4}^{2}\pmod{n}$

implies that

$\displaystyle\sum(x^{8}+6x^{6}z^{2}+6x^{6}y^{2})$ $\displaystyle+$
$\displaystyle\sum(x^{4}y^{4}+x^{4}z^{4}+y^{4}z^{4})$ $\displaystyle+$
$\displaystyle\sum(6x^{2}y^{2}z^{4}+6x^{2}y^{4}z^{2}+36x^{4}y^{2}z^{2})\equiv
S_{4}^{2}\pmod{n}.$

Of these three terms, the first vanishes as discussed, and the third vanishes
by virtue of $\sum x^{4}y^{2}z^{2}\equiv 0\pmod{n}$ and variants. The second
term is $3\cdot S_{4}^{2}\pmod{n}$, whence we arrive at a contradiction,
$3S_{4}^{2}\equiv S_{4}^{2}\pmod{n}$. ∎

However, it is possible to implement our constructions in higher dimensions.

######  Theorem 13.

Let $n,d$ be positive integers. There is a constant $c$, depending on $d$,
such that one can always place $n^{d-1}-cn^{d-2}$ queens on
$\mathbb{Z}_{n}^{d}$ without conflict.

###### Proof.

Clearly, two queens are in conflict if and only if they occupy the same
diagonal. Let us describe a line in $\mathbb{Z}_{n}^{d}$ parallel to
$(\epsilon_{1},\dots,\epsilon_{d})\in\\{-1,0,1\\}^{d}\setminus\\{\mathbf{0}\\}$
as the fiber of a linear map. If we move $\mathbf{x}$ along this diagonal, the
coordinates $x_{i}$ with $\epsilon_{i}=0$ shall be constant, as well as the
sums or differences $x_{j}\pm x_{l}$, where $l$ is the minimal $i$ such that
$\epsilon_{i}=\pm 1$ and $j$ is another such $i$. It follows that our diagonal
is given by the fiber of a diagonal map, a map in the form

$\mathbf{x}\mapsto(x_{1},\dots,x_{l-1},x_{l+1}+e_{l+1}x_{l},\dots,x_{d}+e_{d}x_{l}),\quad
e_{l+1},\dots,e_{d}\in\\{-1,0,1\\}.$

Therefore, queens that occupy a subset of $\mathbb{Z}_{n}^{d}$ are not in
conflict if and only if all maps of this form are injective. To verify this
property, we may compose with an invertible function
$\mathbb{Z}_{n}^{d-1}\to\mathbb{Z}_{n}^{d-1}$, such as that which from one
coordinate subtracts a multiple of another.

Let us place queens on fields parametrised by $\mathbf{t}\in\mathbb{Z}^{d-1}$,
of the form

$\mathbf{F}(\mathbf{t})=(t_{1},-2t_{1}+t_{2},\dots,-2t_{d-2}+t_{d-1},-2t_{d-1}+\alpha(\mathbf{t}))\in\mathbf{Z}_{n}^{d}.$

In order to check independence, let us evaluate a given diagonal map at
$\mathbf{x}=\mathbf{F}(\mathbf{t})$. For $l<i<d$, the component
$x_{i}+e_{i}x_{l}$ is a linear function of $t_{1},\dots,t_{i-1},t_{i}$, where
$t_{i}$ has coefficient $1$. Also, for $i<l$, this is true of the component
$x_{i}$. For the critical parameter $t_{l}$, however, we shall transfer
dependence to the final component through an elimination process.

To first eliminate $t_{d-1}$, we modify
$x_{d}+e_{d}x_{l}=-2t_{d-1}+\alpha+e_{d}x_{l}$ by adding
$2(x_{d-1}+e_{d-1}x_{l})$. The outcome now involves $t_{d-2}$, instead,
eliminated by adding $4(x_{d-2}+e_{d-2}x_{l})$. Recursing, the final component
ends up equal to

$-2^{d-l}t_{l}+\alpha+\left(\sum_{i=0}^{d-l-1}2^{i}e_{d-i}\right)x_{l},$

where the coefficient of $t_{l}$ is given by
$-a=-2^{d-l}+\sum_{i=0}^{d-l-1}2^{i}e_{d-i}$. Our key claim is that for a
suitable step function $\alpha(\mathbf{t})$, for any $0<a<2^{d-l+1}$,
$-at_{l}+\alpha(\mathbf{t})$ is an injective function of $t_{l}$, if a bounded
number of possible values are dismissed.

The task of specifying $\alpha(t_{1},t_{2},\dots,t_{d})$ is accomplished as
follows. Given $l$ and a positive prime $p$, let $q$ be a prime power of $p$
less than $2^{d-l+1}$. If $p$ divides $n$, we subdivide $\\{0,1,\dots,n-1\\}$
at multiples of $\frac{n}{q}$ into $q$ equal intervals, and require that
$\alpha\equiv s\pmod{q}$ for the $s$-th interval in some order. This mesh of
conditions is met by refining modulo $p$, then modulo $p^{2}$, and so on. In
the iteration from $p^{h}$ to $p^{h+1}$, if there are $D$ many coordinates to
refine, we add the chessboard pattern of
$\mathbf{Z}_{p}^{D}\to\mathbf{Z}_{p}$, $(y_{1},\dots,y_{D})\mapsto
y_{1}+\dots+y_{D}$ times $p^{h}$, stretched by $\frac{n}{p^{h+1}}$,
$(p^{h})^{D}$ times. If $p$ does not divide $n$, we simply demand
$\alpha\equiv 0\pmod{q}$. Applying the Chinese Remainder Theorem, the
resulting $\alpha(\mathbf{t})$ is bounded by the product of maximal
$q$.111This product is at most $2^{2^{d+2}}$, by a standard Erdős bound that
considers $\binom{2N}{N}$.

How is $-at_{l}+\alpha$ injective? From the prime powers in the standard form
of $a$, let $c$ be the product of those coprime with $n$, so that
$\frac{\alpha}{c}=\beta$ is integer-valued, yielding
$-at_{l}+\alpha=-c(bt_{l}-\beta)$. If $q$ is a prime power in
$\text{gcd}(b,n)$, $\beta\pmod{q}$ is constant on translates of
$\left\\{0,1,\dots,\frac{n}{q}-1\right\\}$ by multiples of $\frac{n}{q}$. If
$\beta$ was constant, $bt_{l}-\beta$ would attain all of a specific modulo $q$
residue on a suitable translate. A singularity of $\beta$ shifts these values,
still distinct if we drop as many $t_{l}$ values as the jump.

Thus, the range is reduced by $n^{d-2}$ a bounded number of times, resulting
in at least $n^{d-1}-cn^{d-2}$ queens, where $c(d)\leq d\cdot 2^{2^{d+2}}\cdot
2^{2^{d+2}}$. Injectivity follows by successively determining $t_{l}$,
$t_{d-1},\dots,t_{l+1}$, and $t_{l-1},\dots,t_{1}$. ∎

## References

  * [1] Hedetniemi, J. T. and Hedetniemi, S. T. Domination in Chessboards. In Haynes, T. W. et al. (eds.), Structures of Domination in Graphs, Developments in Mathematics 66, 341-386, 2020.
  * [2] Barr, J. and Rao, S. The $n$-Queens Problem in Higher Dimensions. Arxiv preprint, arXiv:0712.2309, 2018.
  * [3] Pólya, H. Über die "doppelt-periodischen" Lösungen des $n$-Damen-Problems. In Ahrens, W. (ed.), Mathematische Unterhaltungen und Spiele 2, B. G. Teubner, 364-374, 1918.
  * [4] Burger, A. P., Cockayne, E. J., Mynhardt, C. M. Queens Graphs for Chessboards on the Torus. Australasian Journal of Combinatorics 24, 231-246, 2001.

Department of Pure Mathematics and Mathematical Statistics, University of
Cambridge, Wilberforce Road, Cambridge CB3 0WB.

E-mail address<EMAIL_ADDRESS>


# A two-level machine learning framework for predictive maintenance:
comparison of learning formulations

Valentin Hamaide  Denis Joassin  Lauriane Castin  François Glineur  ICTEAM
Institute, UCLouvain, Louvain-la-Neuve, Belgium  Ion Beam Applications,
Louvain-la-Neuve, Belgium

###### Abstract

Predicting incoming failures and scheduling maintenance based on sensors
information in industrial machines is increasingly important to avoid downtime
and machine failure. Different machine learning formulations can be used to
solve the predictive maintenance problem. However, many of the approaches
studied in the literature are not directly applicable to real-life scenarios.
Indeed, many of those approaches usually either rely on labelled machine
malfunctions in the case of classification and fault detection, or rely on
finding a monotonic health indicator on which a prediction can be made in the
case of regression and remaining useful life estimation, which is not always
feasible. Moreover, the decision-making part of the problem is not always
studied in conjunction with the prediction phase. This paper aims to design
and compare different formulations for predictive maintenance in a two-level
framework and design metrics that quantify both the failure detection
performance as well as the timing of the maintenance decision. The first level
is responsible for building a health indicator by aggregating features using a
learning algorithm. The second level consists of a decision-making system that
can trigger an alarm based on this health indicator. Three degrees of
refinements are compared in the first level of the framework, from simple
threshold-based univariate predictive technique to supervised learning methods
based on the remaining time before failure. We choose to use the Support
Vector Machine (SVM) and its variations as the common algorithm used in all
the formulations. We apply and compare the different strategies on a real-
world rotating machine case study and observe that while a simple model can
already perform well, more sophisticated refinements enhance the predictions
for well-chosen parameters.

t1t1footnotetext: The research project was funded by the Walloon region under
the BiDMed project.

##  1 Introduction

Predictive maintenance (PdM), or condition-based maintenance, consists of
recommending maintenance decisions based on the information collected through
condition monitoring, usually in the form of time series. It is usually
formulated in one of the two following ways: i) detecting that the machine
under monitoring has entered a faulty state, and therefore predicting that a
failure is coming, or ii) predicting the remaining useful life (RUL) of the
machine. In the scientific literature, those two approaches are referred to as
_diagnostics_ and _prognostics_ respectively. _Prognostics_ is defined by
Jardine et al. [1] in their review as "to predict faults or failures before
they occur" and _Diagnostics_ as "focusing on detection, isolation and
identification of faults when they occur" or as "a procedure of mapping the
information obtained in the measurement/features space to machine faults in
the fault space", i.e. pattern recognition. However, in this paper, those two
concepts are not differentiated as our purpose is to schedule a maintenance
procedure based on the current monitoring information, regardless of the type
of method used. We propose a framework divided into two levels, the first
level consisting of a machine learning model mapping a set of features into a
health indicator and the second level being responsible for the actual
decision making, where an alarm is raised if the health indicator of the first
level crosses a threshold whose value can be optimized.

Predictive maintenance methods can be categorized into model-based,
statistical and machine learning approaches. In this paper, we will focus on
machine learning (ML) approaches. For surveys on model-based and statistical
methods, the reader can refer to [2, 3].

###  1.1 ML based approaches for predictive maintenance

Machine learning techniques applied to fault diagnosis were recently reviewed
in [4]. Those techniques consist of designing machine learning models to
establish the relationship between selected features and the health state of
machines. In [4], those methods are described chronologically, starting with
traditional ML theories that include data collection, feature extraction and
selection and health state recognition with models such as artificial neural
networks (ANN), support vector machine (SVM) and k-nearest neighbors (KNN). In
the present, deep learning methods such as convolutions neural network (CNN),
auto-encoders (AE) or deep belief networks are used, which are capable of
learning both the task and the feature space simultaneously. A common
assumption of those methods is the availability of sufficient labelled data,
which is usually hard to obtain in real-world engineering scenarios [4]. A
more reasonable assumption would be a case where data has been collected for a
few machines that have gone through a failure (or at least a deteriorate
state) and underwent a corrective maintenance as well as for some other
machines that have been replaced without failure (preventive maintenance). In
that case, although different health states are available in the data, an
exact labelling is not available since one does not necessarily know when a
machine has entered a faulty state. In this paper, we investigate different
formulations on how to express such problems with a machine learning
formulation with variations of the support vector machine (SVM) algorithm. The
reason for choosing this algorithm over deep learning approaches is threefold.
First, we want to isolate the formulation or labelling scenario as much as
possible from the algorithm for our comparison. Second, our case study
involves only a few instances of failures and SVM tends to have good
generalization capabilities even with few instances [5]. Finally, the goal of
this paper is not to seek for the best ML algorithm but to compare learning
formulations.

###  1.2 Classification

In the case of _classification_ , we seek to find whether or not a machine
will go to a failure state within a given time-window. The duration of the
time-window is then a parameter that has to be chosen by the user.
Surprisingly, very few papers treat the case of finding incoming faults in
machines with a labelling purely based on failure time (e.g. the last 5 days
are labelled as faulty). One of the papers that is using this kind of approach
is [6] where the authors use a binary classification algorithm with different
horizons to define the faulty class and select the one which optimize a custom
cost that is a trade-off between fault detection and unexploited life (i.e.
replacing the component too soon and not exploiting the full life of the
component).

###  1.3 Anomaly detection

Besides supervised learning, _anomaly detection_ can be used for fault
diagnosis. In the simple univariate case, a threshold can be learned on a
specific feature which is considered anomalous when the threshold is reached
(e.g. in [7]). In the multivariate case, a Gaussian distribution can be fitted
on healthy data and the Mahalanobis distance used as an indicator of health.
This is the approach taken by [8, 9]. Another option is to use the one-class
SVM algorithm with a model learned on healthy data such as the authors did in
[10, 11]. Deep learning methods based on autoencoders were also applied for
anomaly detection in the context of fault diagnosis in [12].

###  1.4 Regression and remaining useful life estimation

Prognostics, the second type of PdM approach, is usually implemented as a
Remaining Useful Lifetime (RUL) estimation [1], that is predicting at each
time step the number of hours/days/cycles remaining before the machine goes
into failure. Various review papers have been published over the past 15 years
on the subject [1, 13, 14, 15, 16, 17] and most of the papers do not
differentiate prognostics from RUL prediction.

RUL estimation is a regression problem. However, it is traditionally not
implemented as a regular supervised learning problem where a mapping is
learned between the sensor inputs and the RUL. The standard approach is to
extrapolate the trend of one or several health indicators previously extracted
from the data via signal processing or machine learning methods until it
reaches a predefined threshold. A systematic approach to RUL prediction in
that sense is presented in [17]. Examples of machine learning approaches using
this kind of framework can be found in [18]. However, a few papers try to
directly map the inputs to the actual RUL. This is the case in [19] where the
authors apply a support vector regression algorithm to the NASA’s Turbofan
engine degradation dataset (CMAPSS) [20] and more recently with deep learning
such as in [21, 22] on the same dataset. However, this dataset is a simulated
dataset with manually introduced failures [17], which can be quite different
from real-life scenarios. In practice, it is usually hard to estimate the RUL
in engineering scenarios, especially in the absence of clear degradation
trends. In this paper, a mapping is learned between the feature space and the
RUL; however, the final goal is not the RUL itself but the decision taken upon
it, i.e. when an alarm for replacement should be raised based on the predicted
RUL. Other regression approaches using a different labelling scheme are also
tested in the paper and discussed in Section 3.3.3.

###  1.5 Decision making

The use of the ML algorithm in the first level of our framework has a
different purpose from what is commonly done within the predictive maintenance
community. Indeed, we are interested to extract a single continuous health
indicator from the learning algorithm rather than identifying a certain state
in the case of classification, an outlier in the case of anomaly detection or
a RUL in the case of regression. This first level is completely unaware of
time, its main purpose is pattern recognition. The resulting health indicator
is then fed to the second level of the framework which is responsible for
taking a decision for maintenance based on a smart aggregation of the past
values of the health indicator.

The rest of the paper is structured as follows: in Section 2, the framework
for predictive maintenance is outlined. In Section 3, various formulations for
the learning problem are described. The decision-making process is presented
in Section 4. In Section 5, the method of validation of the models is
presented. The different methods are tested on a rotating machine case study
in Section 6 and finally Section 7 concludes.

##  2 Predictive maintenance framework

Figure 1: Predictive maintenance framework

A flowchart for our predictive maintenance framework is depicted in Figure 1.
First, data need to be acquired via sensor measurements in the form of time-
series. The purpose of the data acquisition step is to obtain data from
several machines that got a corrective maintenance, i.e. that went through
failure (or at least until a deteriorated state) but also possibly from
machines that were preventively replaced and did not go into a failure state.
We call the data collected from a particular machine from the installation to
the replacement (with or without failure) a _run_. Let us define
$\mathbf{F}^{(r)}\in\mathbb{R}^{n_{r}\times N}$, the feature matrix of the run
$r$, where $n_{r}$ is the number of samples in the time series and $N$ the
number of features. The dataset consists of $R$ runs, i.e. $r=1,...,R$. The
notation $\mathbf{F}^{(r)}(t)$ is used to access the $t^{th}$ sample in the
matrix. We also define $x=\left[F^{(r)}\right]_{r=1,...,R}$ the merged feature
matrix of size $\sum_{r}n_{r}\times N$, and $x_{i}\in\mathbb{R}^{N}$ refers to
sample number $i$ (regardless of time) of the dataset. Some preprocessing is
then achieved on the data collected, which can include filling missing values,
removing sensor faults, normalization, etc.

The next step consists of extracting meaningful information from the raw
sensor measurements $\mathbf{D}^{(r)}$. This is called feature extraction by
the machine learning community while sometimes referred as health indicator
construction [17] in the predictive maintenance literature. In the context of
time-series data, it consists of designing features aggregating sensor inputs
that summarize a certain time period rather than a point in time. Those
aggregations can be done in the time-domain, frequency-domain or time-
frequency domain. A description of the feature design is shown in Figure 2. In
the case of vibration data, the reader can refer to [23] for a review of the
common features that are used for predictive maintenance.

Figure 2: Feature design

A large number of features can be initially produced by the feature extraction
step. Selecting features targets two goals: removing uninformative features
and reducing the size of the problems to mitigate some possible overfitting
and improve computational efficiency. Two main strategies can be used to
select features: a wrapper approach and a filter approach. The wrapper
approach selects the features according to the model in an iterative loop
where features are added or removed according to the score obtained with the
model. The filter approach is model agnostic and selects the features _a
priori_ according to a certain criterion. Since we are comparing different
models, we will choose a filter approach so that we have a set of features
common in all experiments. However, the criteria to select features for the
filter approach should not be supervised, since different labelling strategies
will be tested and choosing one would unfairly favor the algorithm that uses
this particular labelling. Instead, we use an unsupervised minimum redundancy
maximum relevance feature selection for predictive maintenance which uses a
combination of three prognostic metrics to quantify the relevance of a feature
while minimizing redundancy between features [24]. The prognostic metrics
evaluate what a good prognostic parameter should be, i.e. monotonic with
respect to time, behaving in the same ways for different machines and easily
separable between the starting and failure values. For details on the feature
selection method, the reader can refer to [24].

From the selected set of features, we wish to establish a causal relationship
with the health state of the machine. If the current health state is
considered faulty, an alarm should be triggered for repair or replacement.
This is done in two steps. In a first step, the input features are fed to a
model learned on previous machines, which outputs a single number at each time
step (or a single time-series across time), i.e.
$h^{(r)}(t)=f(\mathbf{F}^{(r)}(t))$. Then, in a second step, a decision is
taken upon this output. An alarm is triggered if $g(h^{(r)}(t))$ exceeds (or
goes below) a threshold whose value was optimized beforehand on different
machines, where $g$ is a certain time-window transformation applied to the
health indicator computed from the model. The function $g$ can be the identity
function or a certain aggregation of the output across time such as a moving
average or exponential moving average. This two-level predictive maintenance
approach is depicted in Figure 3.

The next section compares the different strategies that map the set of
features into a single health indicator, i.e. the ML model part in the diagram
of Figure 1.

(a) First level: Three refinement of the ML model

(b) Second level: Decision maker

Figure 3: Two-level predictive maintenance: (a) a model translates a set of
features into a single indicator (b) a decision is made to trigger an alarm
when the health indicator exceeds an optimized threshold (in this case $g$ is
the identity function).

##  3 First level: Problem formulations

In this section, we compare three different strategies to map a set of
features into a health indicator. A first possibility is to simply follow a
single feature, a second requires to train the algorithm on healthy data and
consider as anomalous what deviates from the norm, and finally the most
refined technique is to train on both healthy and unhealthy data with
supervised learning algorithms, which can be formulated either as binary
classification, multi-class classification or regression. As mentioned
earlier, the machine learning models used is this paper are variation of the
support vector machine algorithm (SVM). For the anomaly detection, the one-
class SVM is used, for the classification and multi-class classification, the
standard SVM is used and for the regression, the support vector regression
(SVR) algorithm is used.

###  3.1 Univariate model

In the univariate case, a single feature is followed across time. The selected
feature can be chosen according to engineering expertise or according to the
maximal relevance score obtained by a feature selection algorithm. In our
case, the relevance criterion is the average of three prognostic metrics:
monotonicity, trendability, prognosability [24]. Additionally, some
aggregations can be performed on a certain time-window up to the current time
instant, such as taking the moving average across several hours or even days
if needed. More complex aggregations can also be performed such as exponential
moving average, autoregressive models, etc. The health indictor computed from
the ML model is then

$h^{(r)}=f_{U}\left(\mathbf{F}^{(r)}\right)$

where $f_{U}$ is the model (any type of aggregation or simply the identity
mapping) where $\mathbf{F}^{(r)}$ contains only one feature in this case.

###  3.2 Multivariate anomaly detection

When performing anomaly detection, the algorithm is only trained on healthy
data. When testing on new data, a sample is either marked as inlier or outlier
(i.e. anomalous). In practice, the decision is not binary but is taken based
on a threshold for a decision function. For instance, in the case of anomaly
detection based on a multivariate gaussian distribution fitted on healthy
data, a decision to flag a sample as anomalous is taken if it is far from the
fitted n-dimensional ellipsoid center. Usually, the distance measure chosen is
the Mahalonibis distance between the sample and the ellipsoid center. The
evolution of that distance across time will therefore be the health indicator
on which a decision to trigger an alarm will be made. In our approach, the
algorithm chosen for anomaly detection is the one-class SVM, rather than an
approach based on the multivariate gaussian algorithm. Even though the
Mahalanobis distance does not apply to the one class SVM, there is a similar
idea of distance from a set of _normal_ samples.

One-class SVM was proposed as an extension of the support vector machine to
the case of unlabeled data [25]. It tries to estimate the distribution of the
input data (considered healthy) by a simpler subset of the input space and
estimates a function $f$ that is positive in that subset and negative on the
complement. The corresponding model is formulated as follows:

$\displaystyle\min_{w,b,\zeta}\frac{1}{2}\lVert
w\rVert^{2}+\frac{1}{\nu\ell}\sum_{i=1}^{n}\zeta_{i}-b$ (1)
$\displaystyle\textrm{s.t. }(w^{T}\phi(x_{i}))\geq b-\zeta_{i},$ (2)
$\displaystyle\zeta_{i}\geq 0,i=1,..,n$ (3)

where $x_{i}$ are the training vectors, $w$ and $b$ the weights and bias for
which we solve, $\zeta_{i}$ are the slack variables allowing some samples to
be on the _wrong_ side of the hypersurface, $\phi(\cdot)$ is a non-linear
mapping to allow for a non-linear boundary, $\ell$ is the number of samples in
the training set and $\nu\in[0,1)$ is a hyperparameter representing an upper
bound on the fraction of training errors and a lower bound of the fraction of
support vectors. For our experiments, as well as the other SVM-based models,
we used the sci-kit learn implementation which is wrapper around the LIBSVM
library [26]. The decision-making process is based on the distance to the
hypersurface function, and the decision function is defined as

$f_{\mathrm{1SVM}}(x)=\sum_{i=1}^{n}\alpha_{i}K(x_{i},x)+\rho$ (4)

where $K(x_{i},x)=\phi(x_{i})^{T}\phi(x)$ is the kernel function and
$\alpha_{i}$ and $\rho$ the dual variable and independent term of the
optimization problem (1). The health indicator computed from the model is then

$h^{(r)}=f_{\mathrm{1SVM}}\left(\mathbf{F}^{(r)}\right)$

The same notation will be used throughout the paper for all SVM-based models
for variables, weights, slacks and kernels.

###  3.3 Supervised learning

Supervised learning refers to algorithms that learn a function mapping from a
set of input variables, the features, to a corresponding output, the labels.
This requires having the corresponding desired label for each input sample.
Those labels can either be numeric values in case of regression or categorical
variables in case of classification. In the context of predictive maintenance,
those categorical variables can be the fact that a machine is in a healthy or
unhealthy state or a certain type of faults on the machine. If the labels are
known, for instance with a healthy machine and a faulty machine on a test
bench, then the problem becomes simple and supervised learning is the way to
go. However, labels are usually not available in real-case scenarios. Indeed,
we do not necessarily know when exactly a machine enters a faulty state, even
if the machine goes into failure at the end of the run. Instead, we use a
labelling based on the remaining time before failure. In the case of
classification and multi-class classifications, the labels are chosen somewhat
arbitrarily by splitting the run into two or more time-periods representing
healthy or unhealthy states. The chosen duration used to split the run into
classes is a parameter that needs to be chosen by the user based on
engineering expertise or tuned as an hyperparameter of the problem. In the
case of regression, the time before failure can be directly used as the
labelling but other possibilities exist and are described in Section 3.3.3.

####  3.3.1 Binary classification

In the classification approach, a run is divided into two classes: faulty (F)
and non-faulty (NF) as follows:

$y^{(r)}(t)=\begin{cases}NF&\text{if $t<T^{(r)}-w$}\\\ F&\text{if $t\geq
T^{(r)}-w$}\end{cases}$

where $T^{(r)}$ is the duration of the run $r$ and $w$ is a parameter to
choose for the faulty state duration. We thus have to choose a duration _a
priori_ taking into account that a $w$ too small will lead to a late detection
while a $w$ too big could lead to too early detection (and therefore
unexploited lifetime). In the case of a machine preventively replaced (no
actual failure at the end of the run), the entire run is marked as $NF$.

The classification algorithm used here is the well-known support vector
machine algorithm [27] which finds the hypersurface that separates the classes
with maximal margin by solving the following optimization problem:

$\displaystyle\min_{w,b,\zeta}\frac{1}{2}\lVert
w\rVert^{2}+C\sum_{i=1}^{n}\zeta_{i}$ (5) $\displaystyle\textrm{s.t.
}y_{i}(w^{T}\phi(x_{i})+b)\geq 1-\zeta_{i},$ (6) $\displaystyle\zeta_{i}\geq
0,i=1,..,n$ (7)

where $y_{i}\in\\{-1,1\\}^{n}$ are the target values (where $F$ is mapped to 1
and $NF$ to -1) and $C$ is an hyperparameter representing the trade-off
between the margin width and the sum of the slack variables
($\sum_{i}\zeta_{i}$).

The decision function is based on the distance to the hypersurface, that is
defined as

$f_{\mathrm{SVM}}(x)=\sum_{i=1}^{n}y_{i}\alpha_{i}K(x_{i},x)+\rho$ (8)

and the computed health indicator is thus

$h^{(r)}=f_{\mathrm{SVM}}\left(\mathbf{F}^{(r)}\right)$

####  3.3.2 Multi-class classification

In the multi-class classification approach, a run is separated into more than
two classes. Each class represents a different non-overlapping time window
between the beginning of the run and the failure. Mathematically, the labels
are defined as

$y^{(r)}(t)=\begin{cases}NF&\text{if $t<T^{(r)}-w_{1}$}\\\ F_{1}&\text{if
$T^{(r)}-w_{1}\leq t<T^{(r)}-w_{2}$}\\\ F_{2}&\text{if $T^{(r)}-w_{2}\leq
t<T^{(r)}-w_{3}$}\\\ \vdots\\\ F_{N-1}&\text{if $t\geq
T^{(r)}-w_{N}$}\end{cases}$

where $NF$ is considered the healthy class and $F_{i},i=1,...,N-1$ are
considered the $N-1$ faulty classes with $w_{1}>w_{2}>...>w_{N}$ and increased
level of fault severity. In the case of a machine preventively replaced (no
failure at the end), the entire run is marked as $NF$.

The multi-class classification approach is similar to the binary approach. The
SVM algorithm is still used but instead of solving one optimization problem as
in (5-7), we solve $N$ optimization problems where $N$ is the number of
classes. We use a one-versus-rest strategy where we train a single classifier
per class, with the samples of that class being positive ($+1$) and all the
rest being negative ($-1$) and thus keep the same formulation as in (5).

Since there are $N$ optimization problems, there are $N$ decision functions
such as in (8). The decision functions are defined as

$h_{j}^{(r)}=f_{\mathrm{SVM}}^{j}\left(\mathbf{F}^{(r)}\right)\quad\text{for
}j=1,...,N$

####  3.3.3 Regression

In the regression approach, three labelling scenarios are tested. The first
approach is to directly use the remaining useful life (RUL) as labels. In the
second approach, instead of using absolute times, relative times are used with
the percentage of life used as labels. Finally, a third approach tries to
mimic the intuition that a machine is stable in the beginning of its life and
deteriorates more and more starting some time before the failure with a
piecewise linear function. We refer to this approach as _ReLu_ , an analogy to
the rectified linear unit in machine learning due to the shape of the
labelling function. The three labelling approaches are detailed in Table 1.
Those labelling strategies are only valid for corrective maintenance. Indeed,
runs of machines preventively replaced have to be disregarded in case of RUL
or percentage of life strategies. For the ReLu strategy, runs of machines
preventively replaced can be kept and labelled as zero across the entire life
span of those machines.

|  RUL |  Percentage of life |  ReLu
---|---|---|---
Labelling |  $h^{(r)}(t)=\frac{T^{(r)}-t}{D}$ with $D$ a normalizing constant |  $h^{(r)}(t)=\frac{t-t_{0}^{(r)}}{T^{(r)}-t_{0}^{(r)}}$ where $t_{0}$ is the start time and $T$ the end time. |  $\displaystyle h^{(r)}(t)=$ $\displaystyle\begin{cases}0&\text{if $t<T-t_{d}$}\\\ \frac{T-t}{t_{d}}&\text{if $t\geq T-t_{d}$}\end{cases}$ where $T$ the end time and $t_{d}$ is a supposed start of deterioration fixed for training.
Function shape |  T0 |  01 |  01
Decision |  $g\left(h^{(r)}\right)<L$ |  $g\left(h^{(r)}\right)>L$ |  $g\left(h^{(r)}\right)>L$
Table 1: Labelling strategies for regression & decision making

Support vector regression (SVR) is the common regression algorithm used for
all labelling scenarios. It is a variation of the SVM performing a regression
based on the concept of support vectors [28]. The idea is to find a function
$f(x)$ that has at most an $\varepsilon$ deviation from the targets $y_{i}$ by
solving the following optimization problem

$\displaystyle\min_{w,b,\zeta,\zeta^{*}}\frac{1}{2}\lVert
w\rVert^{2}+C\sum_{i=1}^{n}(\zeta_{i}+\zeta_{i}^{*})$ (9)
$\displaystyle\textrm{s.t.
}y_{i}-w^{T}\phi(x_{i})-b\leq\varepsilon+\zeta_{i},$ (10) $\displaystyle
w^{T}\phi(x_{i})+b-y_{i}\leq\varepsilon+\zeta_{i}^{*},$ (11)
$\displaystyle\zeta_{i},\zeta_{i}^{*}\geq 0,i=1,..,n\ $ (12)

where $y_{i}$ are the targets values, $\zeta_{i}$ and $\zeta_{i}^{*}$ are the
slack variables allowing some samples to be outside of the tube of radius
$\varepsilon$ centered around the function and $C$ is an hyperparameter
representing the trade-off between the flatness of $f$ and the sum of
deviations larger than $\varepsilon$ (that is
$\sum_{i}\zeta_{i}$+$\zeta_{i}^{*}$). For more information on the SVR and how
this optimization problem can be solved efficiently, the reader can refer to
[29].

The decision function of the SVR is defined as

$f_{\mathrm{SVR}}(x)=\sum_{i=1}^{n}(\alpha_{i}-\alpha_{i}^{*})K(x_{i},x)+\rho$
(13)

where $\alpha_{i}$, $\alpha_{i}^{*}$ and $\rho$ are the dual variables and
independent term of the optimization problem (9). The computed health
indicator (HI) is thus

$h^{(r)}=f_{\mathrm{SVR}}\left(\mathbf{F}^{(r)}\right)$

##  4 Second level: Decision making

The purpose of the second level of the predictive maintenance framework is the
make a decision on the health indicator computed from the ML model. Let
$h^{(r)}=f(\mathbf{F}^{(r)})$ be the HI computed from the model, the decision
maker will trigger an alarm for replacement if $z=g(h^{(r)})>L$ where $g$ is a
function of time-series $h^{(r)}$ and $L$ is a predefined threshold. In the
simplest form, $g$ is the identity function, or the negative function
$g(h^{(r)})=-h^{(r)}$ in case $L$ is a lower bound. When the current output is
above (or below) the threshold $L$, an alarm is raised. The function $g$ can
also be an aggregation of past values such as a moving average or an
exponential moving average. In case $g$ is a moving average, it is defined as

$g(y)=\frac{h(t)+h(t-1)+...+h(t-H)}{H}$ (14)

where $H$ is the horizon selected for the aggregation. In case $g$ is an
exponential moving average, the aggregation is defined recursively as

$\displaystyle z_{0}$ $\displaystyle=$ $\displaystyle h_{0}$ (15)
$\displaystyle z_{t}$ $\displaystyle=$ $\displaystyle\alpha
h_{t}+(1-\alpha)h_{t-1}$ (16)

The difference with the simple moving average is that in this case, the
window-size is infinite but the weights are exponentially decreasing. However,
the $\alpha$ parameter can be tuned to be interpreted approximately as a
H-hour moving average when computed as

$\alpha=\frac{2}{H+1}$ (17)

where $H$ is the horizon111This is an approximation of of the formula
$\alpha=1-\exp(\log(1-p)/H)$ where $p=0.86$ is the contribution of the window
on the moving average. For instance, if $H=12$ and the time between two
consecutive sample is 1 hour, $\alpha=\frac{2}{12+1}\approx 0.1538$ and is
interpreted as a 12-hour moving average.

In the case of binary classification, an alarm would be triggered if $z(t)>L$
where $L$ is a previously optimized threshold. For regression, the idea is
similar and the decision functions are detailed in Table 1. For the multi-
class classification, since there are multiple decision functions, the process
has to be slightly adapted. An alarm is triggered if the value of the decision
function of one of the faulty classes ($F_{1},...,F_{N-1}$) is higher than the
healthy class $NF$. Mathematically, we trigger an alarm at sample $i$ if any
$j\neq 1$ satisfies

$g\left(f_{\mathrm{SVM}}^{F_{j}}\left(\mathbf{F}^{(r)}\right)\right)>g\left(f_{\mathrm{SVM}}^{NF}\left(\mathbf{F}^{(r)}\right)\right)$
(18)

where $F_{j}$ is the decision function associated to class $F_{j}$,
$j=1,...,N-1$, $NF$ the decision function associated to the healthy class and
$g$ is the function applied to the health indicator.

##  5 Assessing the predictive performance of models

###  5.1 Scoring

The end goal of a predictive maintenance application is to help on the
decision to replace or repair the machine under monitoring at the right time.
The right time may vary between applications but is usually a trade-off
between detecting failure and limiting false alarms. As a first goal, we want
to check whether or not a fault (in the case of a run with failure) can be
detected and minimize false alarm. As a second goal, we want to score the
timing of the alarm. This concept can be translated into two metrics. We
define a false positive as an alarm that was triggered too early, in our case
more than 15 days in advance. A true positive is defined as an alarm raised
between 0 and 15 days in advance. While true positives are only relevant for
corrective maintenance, false positives are relevant for both corrective and
preventive maintenance. Indeed, an alarm raised at any time for a machine
which was preventively replaced is considered a false positive. We thus define
two metrics, the false positive rate:

$\text{FPR}=\frac{\text{FP}}{C+P}$ (19)

and the true positive rate

$\text{TPR}=\frac{\text{TP}}{C}$ (20)

where $C$ and $P$ are the numbers of corrective and preventive runs
respectively, and FP and TP the numbers of false positives and true positives
respectively. Note that we do not assess each individual prediction for each
time step, but only evaluate the quality of the first trigger for each run
(hence the earliest for that run).

We combine those two metrics into a single metric by taking their harmonic
average, similarly that we would do to compute the $F_{1}$ score between
precision and recall in a conventional classification scenario. Moreover, we
add a $\beta$ parameter, controlling the importance of the TPR over the
$1-\text{FPR}$. We call this metric, the $F_{\text{score}}$ and define it as

$F_{\text{score}}=(1+\beta^{2})\cdot\frac{(1-\text{FPR)}\cdot\mathrm{TPR}}{\beta^{2}\cdot(1-\mathrm{FPR})+\mathrm{TPR}}$
(21)

Since it is also important to score the quality of timing at which the alarm
is raised, another metric called the business metric ($B_{\text{score}}$) is
defined. It provides a score between zero and one (higher is better) computed
from a piecewise-linear function of the number of days between the prediction
and the actual failure:

  * 1.

A prediction 7 days ahead is considered optimal and gives a perfect score.

  * 2.

Predictions between 7 and 0 days assign a score that decreases linearly to
zero (with a slightly lower slope between 7 and 5 days).

  * 3.

Predictions ranging from 7 to 15 days are assigned a score linearly decreasing
from 1 to 0 (with a slightly lower slope between 7 and 10 days).

  * 4.

Predicting a failure more than 15 days in advance leads to a zero score to
reflect the unexploited lifetime.

Note that this business score is only applicable for corrective maintenance.
While $F_{\text{score}}$ is a single score resulting from an ensemble of runs,
$B_{\text{score}}$ is defined per run and one must take the average across all
corrective maintenance to obtain a single score, i.e.
$\overline{B_{\text{score}}}=\sum_{r}\frac{B_{\text{score}}^{(r)}}{C}$.
Finally, we combine the $F_{\text{score}}$ and the business score into a
single score that takes into account both the corrective and preventive
maintenance runs. We define $\alpha\in[0,1]$ as the weight associated to the
$F_{\text{score}}$ and $1-\alpha$ the weight associated to
$\overline{B_{\text{score}}}$. Since the $F_{\text{score}}$ is used for both
corrective and preventive maintenance while $B_{\text{score}}$ is only applied
to corrective maintenance, we must further multiply $1-\alpha$ (the weight
associated to $B_{\text{score}}$) by the ratio of corrective maintenance over
the the number of runs, i.e. $\frac{C}{C+P}$. In the end, the final score is
defined as

$\frac{\alpha}{(1-\alpha)\frac{C}{C+P}+\alpha}F_{\text{score}}+\frac{(1-\alpha)\frac{C}{C+P}}{(1-\alpha)\frac{C}{C+P}+\alpha}\overline{B_{\text{score}}}$
(22)

###  5.2 Cross validation

In order to obtain an unbiased estimation of the performance of the algorithms
as well as tuning the different hyperparameters of the machine learning model
and the threshold for the decision making, one must cross-validate the
results. The usual way to validate a model in the context of machine learning
is to split the data into a training, validation and test set. The training
set is used to train the algorithm, the validation set is used to tune the
hyperparameters of the model and the test set, a completely independent set,
is used to assess the final prediction on unseen data.

In the context of predictive maintenance, some precautions are necessary.
Splitting the data into training, validation and test sets cannot be done in a
completely random way. Recall that prediction occurs in a continuous fashion
along the time series, i.e. we classify or predict at each time step. Hence,
due to the temporal nature of the prediction task, data are correlated in time
and one cannot use information learned in the future to predict the past or
the present. Therefore, the training set should not contain data that are
further in time than the validation and test set. An even better practice is
to split data per run, meaning that data from a particular machine cannot be
split among different sets. This is the approach that we take. Some machines
are assigned to the training set, some to the validation set and the rest to
the test set.

However, in most real-life applications, machine runs-to-failure are scarce.
Hence, it is difficult to build a sufficiently large (in terms of number of
corrective maintenance) training set or test set. Therefore, we perform a
cross-validation, where a fold is defined as a run. However, a single layer of
cross-validation, for example a leave-one-out cross-validation, is still
biased because the hyperparameter optimization has _seen_ all data. Since data
are scarce, it is not an option to leave a few runs as test set as they would
most likely not represent very well the behaviors of all machines. Instead, we
perform what we call a double cross validation. The framework for double
cross-validation is outlined in Figure 4.

Figure 4: Example of double cross validation process for 9 machines subdivided
into 3 folds.

This double cross-validation consists of two levels of cross-validations. In
the inner loop a leave-one-out cross validation is performed to tune the
hyperparameters of the ML model as well as tune the threshold for the
decision-making. Then, a model is trained on all data of the inner loop
(training set + validation set) with the best set of hyperparameter found in
that inner loop.

In the outer loop, the model is tested on an independent set of runs. Then, a
new test set is selected at the next outer loop. In the end, we obtain not one
but multiple models (with possibly different set of hyperparameters) that
correspond to the number of folds on the outer loop. The average of those
scores results in an unbiased estimation of the performance of our algorithms.
Note that a fold can contains both corrective or preventive maintenance. We
target on an equal distribution of those two categories of runs across the
different folds.

##  6 Application to a rotating machine

###  6.1 Problem description

The predictive maintenance case study we consider in this work deals with a
high-speed rotating condenser (RotCo) modulating the RF frequency inside a
synchrocyclotron [30]. The RotCo is composed of a stator and a rotor with
eight blades and rotates at a constant speed of 7500 RPM with the help of ball
bearings. A picture of the system is shown in Figure 5.

Figure 5: (bottom) RF system of the cyclotron with the rotating condenser on
the right. (Top) Detailed view on the rotating condenser. [30] CC-BY-3.0

Several sensors located inside the machine are used to gather data. An
accelerometer sensor placed on the condenser external casing measures
vibrations and performs 10-second acquisitions at a rate of 10kHz, once every
hour. Four other sensors are placed on the machine to gather data every
second. Those include two temperature sensors, a vacuum pressure and a torque
sensor. In total, 11 corrective maintenance and 28 preventive maintenance runs
have been gathered, i.e. a total of $R=39$ runs.

After this data acquisition step, several features are built. For the
vibration data, time-domain and frequency-domain features on each of the
10-second long acquisition files are constructed. Table 2 details the way
those features have been computed.

Raw signals |  Feature extraction |  Selected features
---|---|---
_Vibration amplitude_ |  _Time-domain:_
RMS, MAD, Peak to Peak Amplitude, Skewness, Kurtosis, Crest Factor, Clearance
Factor, Shape Factor, Margin Factor, Max Amplitude

_Frequency-domain:_
Amplitude 1N2221N refers to one time the fundamental frequency., 2N, 3N, every 20Hz band from 0-1kHz, BPFO 1-3N, BPFI 1-3N, BSF 1-3N, FTF 1-3N |  MAD, Margin Factor, BPFO 3N, Peak to Peak Amplitude, Crest Factor, Amplitude 1N, RMS, Spectral Amplitude at 350 +- 10 Hz
_Non-vibration features_ : Bearing Temperature
Pyrometer Temperature
Torque
Vacuum Pressure |  _Time-domain:_
Mean, Max, Min, Standard Deviation, Skewness, Kurtosis |  Torque Mean, Torque Max
Table 2: Feature design

In the time-domain, those include Root Mean Square (RMS), Median absolute
deviation (MAD) which is a robust measure of variability based on the
deviations from the median, peak to peak values, skewness (third statistical
moment) and kurtosis (fourth statistical moment). Those also include metrics
based on the peaks of the signals: crest factor, clearance factor, shape
factor, margin factor and max amplitude (for a detailed explanation on those
features, the reader can refer to [31]). For the frequency-domain features, we
compute the amplitudes at the fundamental frequency and its first three
harmonics, the spectral power of all 20 Hz non-overlapping bands from 0-1kHz
and finally the amplitudes at characteristic bearing frequencies [32], i.e.

  * 1.

Ball Pass Frequency Outer Race (BPFO):
$\frac{n_{f}}{2}\left(1-\frac{D_{b}}{D_{p}}\cos\phi\right)$

  * 2.

Ball Pass Frequency Inner Race (BPFI):
$\frac{n_{f}}{2}\left(1+\frac{D_{b}}{D_{p}}\cos\phi\right)$

  * 3.

Ball spin frequency (BSF):
$\frac{D_{p}f}{2D_{b}}\left(1-(\frac{D_{b}}{D_{p}}\cos\phi)^{2}\right)$

  * 4.

Fundamental train frequency (FTF):
$\frac{f}{2}\left(1-\frac{D_{b}}{D_{p}}\cos\phi\right)$

where $D_{p}$ is the pitch diameter, $D_{b}$ the ball diameter, $\phi$ the
contact angle and $n$ the number of balls. The first three harmonics of those
characteristic frequencies are also included. For the non-vibration data,
aggregations of the signals over a 1-hour time-window are performed to match
with the vibration acquisition sampling. Those aggregations include the mean,
max, min, standard deviation, skewness and kurtosis values. This finally
results in $N=89$ features (65 from vibration data and 24 from non-vibration
data) computed every hour. The features are then scaled by subtracting their
mean and dividing by their standard deviation, and then smoothed by a moving
average over a 12-hour time window.

The ten best features are then selected according to the feature selection
method presented in [24]. This is an unsupervised feature selection;
therefore, no formulation is favored.

All the formulations presented in Section 3 are tested on this real-world
application. To avoid comparing all combinations of the approach of the first
level with the ones of the second level, we first compare the approach of the
first level with the simplest second level decision-making process, where the
function $g$ is the identity function. This means that the health indicator is
directly compared to a threshold and an alarm for replacement is raised
whenever this output exceeds (or goes below) the optimized threshold. This
allows us to select a subset of formulations obtaining the best scores, which
are then tested against different aggregations $g$ in Section 6.3.

The double cross-validation is performed as follows: 11 corrective maintenance
runs are distributed among 11 folds and the 28 machines preventively replaced
are distributed equally among those 11 folds. After leaving one fold aside as
test set, we perform the inner loop of the cross-validation where the ML model
is trained consecutively on all but one run (where a fold is defined as a
run). The inner cross-validation allows to tune the hyperparameters of the
model. Then the test set changes to the next outer fold and the whole process
starts again.

###  6.2 First level: training & results

Table 3 summarizes the different hyperparameters of the problem for each
formulation as well as the different labelling scenarios tested. Note that the
univariate model is absent from the table because no training is required at
the first level of the framework.

Formulation |  Hyperparameters |  Scoring metric
---|---|---
Binary classification |  $C=[10^{-2},10^{-1},1,10]$
$\gamma=[10^{-5},10^{-4},10^{-3}]$
kernel: linear, RBF
horizon = [3,5,7,10] days |  F-score333The F-score is defined as the harmonic mean between precision and recall, i.e. $F_{\text{score}}=\frac{2\times\text{precision}\times\text{recall}}{\text{precision}+\text{recall}}$
Multi-class classification |  $C=[10^{-2},10^{-1},1,10]$
$\gamma=[10^{-5},10^{-4},10^{-3}]$
kernel: linear, RBF
labelling 1: 3 classes (0-5 days, 5-10days, >10days)
labelling 2: 6 classes (0-2 days, 2-4days, …, 8-10days, >10days) |  Unweighted mean of F-score associated to each label
One-class SVM |  $\nu=[0.01,0.05,0.1,0.5]$
$\gamma=[10^{-4},10^{-3},10^{-2}]$
kernel: linear, RBF
Horizon = 15 days |  F-score (with failure data included during testing)
Regression (RUL, RUL percentage, ReLu) |  $C=[10^{-2},10^{-1},1,10]$
$\gamma=[10^{-5},10^{-4},10^{-3}]$
$\epsilon=[0.01,0.1,0.5]$
kernel: linear, RBF
horizon for ReLu: $t_{d}=10$ days |  Mean absolute error: $\text{MAE}(y,\hat{y})=\frac{1}{n}\sum_{i=0}^{n_{-}1}\left|y_{i}-\hat{y}_{i}\right|.$
Table 3: Training & hyperparameter tuning for first level of the predictive
maintenance framework

For binary classification, several horizons are tested to split the data among
healthy and unhealthy, ranging from 3 to 10 days. The F-score is used to
select the best set of hyperparameters. For multi-class classification, two
labelling scenarios are tested. The first one includes three classes defined
in the following way: from 0 to 5 days prior to failure, 5 to 10 days, and
more than 10 days. The second labelling scenario includes 6 classes defined as
2-day periods from the failure and a class defined as more than 10 days prior
to failure. Hyperparameters are selected according to the mean of the F-score.

For the One-Class SVM implementing anomaly detection, the model is only
trained on data further than 15 days prior to failure for corrective
maintenance and all data for preventive maintenance. To select the best
hyperparameters set however, the model is tested on all data of a run and the
hyperparameters that obtained the best F-score are selected.

For the RUL and RUL percentage formulations, the training can only be done on
corrective maintenance runs. For the ReLu formulation, preventive runs can be
included, as the labels associated to those runs can be labelled as zero
across their lifetime, since the ReLu is defined as nonzero only at a time
$t_{d}$ prior to failure and monotically increasing until the actual failure.
Parameter $t_{d}$ is fixed at 10 days in all our tests. The metric used to
select the hyperparameters is the Mean absolute error for all regression
formulations.

The results of the double cross-validation for all formulations is outlined in
Table 4. The count of true positives and false positives are summed across the
different folds of the test set (every run is at least in one test set) and
the $F_{\text{score}}$ of equation (21) is computed with those counts with
parameter $\beta=0.5$ to give more emphasis on avoiding false positives. The
business score is averaged on all corrective maintenance runs since the
business score is not defined for runs with preventive replacement. Finally,
the final score is computed according to equation (22) with parameter
$\alpha=0.75$ to increase the emphasis on the $F_{\text{score}}$.

Formulation |  False Positives |  True Positives |  Business score |  F score |  Final score
---|---|---|---|---|---
Univariate model |  2/39 |  8/11 |  0.436 |  0.894 |  0.855
One-class SVM |  16/39 |  5/11 |  0.257 |  0.556 |  0.531
Binary classification |  3 days: 6/39
5 days: 5/39
7 days: 4/39
10 days: 2/39 |  3 days: 9/11
5 days: 9/11
7 days: 9/11
10 days: 9/11 |  3 days: 0.495
5 days: 0.526
7 days: 0.465
10 days: 0.500 |  3 days: 0.840
5 days: 0.860
7 days: 0.880
10 days: 0.919 |  3 days: 0.811
5 days: 0.831
7 days: 0.845
10 days: 0.883
Multi-class classification |  3 classes: 5/39
6 classes: 5/39 |  3 classes: 9/11
6 classes: 8/11 |  3 classes: 0.493
6 classes: 0.468 |  3 classes: 0.860
6 classes: 0.838 |  3 classes: 0.829
6 classes: 0.807
RUL |  32/39 |  2/11 |  0.141 |  0.180 |  0.177
RUL percentage |  12/39 |  6/11 |  0.254 |  0.657 |  0.622
ReLu |  6/39 |  10/11 |  0.617 |  0.858 |  0.837
Table 4: Results of first level. A run is considered a false positive if an
alarm is raised more than 15 days in advance. An alarm is considered a true
positive if it was raised between 0 and 15 days prior to failure in case of a
run with failure. The business score is the mean of the business scores for
all corrective maintenance. Bold values are the best results for each
criterion. Final score is computed via equation (22) with $\alpha=0.75$ and
$\beta=0.5$ and is the score of interest that determines the best formulation.

Surprisingly, we observe that the simplest model, the univariate model based
on the best feature with respect to the three prognostic metrics from [24],
performs quite well, better than most of the approaches tested even though the
other approaches also include this feature within their ten selected features.
This feature is the median absolute deviation of the vibration amplitude
averaged over the last 12 hours. It is the median of the absolute deviation
from the data’s median and is computed as
$\text{MAD}(x)=\text{median}(|x_{i}-\bar{x}|)$ where
$\bar{x}=\text{median}(x)$.

In the second degree of refinement (see Figure 3(a)), the one-class SVM
performs poorly with many false positives resulting in a low final score. For
the third refinement, i.e. supervised learning, we can make several
observations. In the binary classification, the more we increase the size of
the faulty class, the better performance we obtain. This could be explained by
two different reasons. The first one is that the class imbalance is reduced
when the size of the faulty class increases. The second reason is that signs
of faults already appear up to 10 days in advance. In the multi-class
classification, performances are slightly lower than for the binary
classification. Thus, splitting the runs into multiple classes does not seem
to help.

For the three regression formulations, the results are quite different from
each other. Directly mapping the features input to the RUL does not seem to
work at all. This could be explained by the fact that there is too much
disparity between the life spans of the different runs, or simply not enough
failure samples. This can also partly be explained by the fact that the
$\epsilon$ parameter in the SVR formulation should be tuned more carefully for
this type of formulation, since the labels are not scaled to 0-1 like for the
RUL percentage or ReLu. The percentage of life formulation performs better
than the conventional RUL but is still far behind the other formulations. This
could also be explained by the fact that training is only performed on
corrective maintenance for those two formulations. Finally, the ReLu
formulation performs quite well although the number of false positive is high.

In conclusion for this analysis of first-level formulations, we find that no
method clearly outperforms all the others and that depending on which criteria
we focus on (i.e. which column of Table 4 we look at), several formulations
can be recommended. The trade-off between high failure detection rate and low
false alarm is one of those determining aspects. In our case study the binary
classification formulation gives the best results in terms of the final score
for a well-chosen window size.

In the next section, different decision-making function $g$ are tested on the
best algorithms obtained at the first level.

###  6.3 Second level: training & results

In this section, we compare the application of a function $g$ on the computed
health indicators of the first level that lead to the best scores, i.e.
classification with a 10 days window for the faulty class, the multi
classification with 3 classes, the ReLu formulation and the univariate model.
We compare the case where $g$ is a moving average and an exponential
smoothing. For all cases, we test a moving window of size 12-hours, 24-hours,
48-hours and 5-days. In the case of the exponential smoothing, since we have
an infinite window, we use equation (17) to match the parameter $\alpha$ with
the window size. The results are outlined in Figure 6 and a more detailed
version is available in the appendix Tables 5-8. We observe that aggregating
the computed HI with a moving average or exponential smoothing is not a
guarantee for better results. However, by tuning the window size parameter, we
are able to obtain better performance than the identity mapping for the
exponential smoothing in all models.

We observe that when the window size at the second level increases, the number
of false alarms (false positives) decreases, but so as the number of cases
detected (true positives) which sometimes translates into a lower score. The
absolute minimization of false alarms might be wanted in some applications;
thus, this second level of aggregation could be effective in that case. We
also observe that the exponential smoothing performs better in general than
the simple moving average. Therefore, using a scheme of decreasing weights
with respect to time might be a good idea. The relatively modest impact of the
function $g$ on the final results might also be due to the fact that the
features inputs were also averaged across a 12-hour time-window before the
model is applied.

In the end, it is hard to give a definitive conclusion about the best
formulation to use in a predictive maintenance scenario, due to the scarcity
of runs and failure data, which greatly impacts the final score. However, some
insights can still be taken. A simple univariate model can already be
effective provided that a health indicator with high predictive power has been
found. If we have corrective maintenance data, a supervised learning algorithm
should lead to better performance than a one-class classification or anomaly
detection algorithm. The binary classification formulation gives the best
results when the time-window splitting healthy and faulty data is carefully
selected. Finally, while directly mapping the RUL to the inputs is not a good
idea, we find that a formulation such as the ReLu formulation which can be
trained on both preventive replacements runs and runs-to-failure with a
labelling mimicking an increasing fault severity starting at a certain time
$t_{d}$ before the failure, results in good performance.

Figure 6: Results of second level

##  7 Conclusion

The aim of predictive maintenance is to avoid failure by replacing or
repairing a machine at the right time. This is done by attempting to raise an
alarm before the failure, and ideally sufficiently in advance to ease the
maintenance scheduling. The best timing for detecting a failure was encoded
via a business metric in our work. We developed a two-level framework to
tackle this problem. In the first level, we compared several formulations to
discriminate healthy and unhealthy states of the machines based on the
remaining time before failure, in the absence of true labelled machine
malfunctions. In the second level, we compared different ways to exploit the
health indicator computed from the learning algorithm and we optimized a
threshold for making the decision on when to raise an alarm. Our two-level
framework approach gives promising results on the rotating machine case-study
we considered. One of the key take-away is that the most complex models do not
necessarily give the best results, and an univariate model can already perform
very well when a powerful predictive feature can be first extracted. Although
more complex multivariate models rely on more information and should
theoretically lead to a more robust solution, it is not always the case.
Special care must be taken to select the proper hyperparameters. Depending on
the sought trade-off between false alarm and failure detection rate, certain
methods and thresholds will perform better than other combinations, which
should be carefully selected by the user and probably depend on the considered
case study.

## References

  * [1] A. K. Jardine, D. Lin, D. Banjevic, A review on machinery diagnostics and prognostics implementing condition-based maintenance, Mechanical systems and signal processing 20 (7) (2006) 1483–1510.
  * [2] Z. Gao, C. Cecati, S. X. Ding, A survey of fault diagnosis and fault-tolerant techniques—part 1 : Fault diagnosis with model-based and signal-based approaches, IEEE transactions on industrial electronics 62 (6) (2015) 3757–3767.
  * [3] Z. Gao, C. Cecati, S. X. Ding, A survey of fault diagnosis and fault-tolerant techniques—part 2 : Fault diagnosis with knowledge-based and hybrid/active approaches, IEEE Transactions on Industrial Electronics 62 (6) (2015) 3768–3774.
  * [4] Y. Lei, B. Yang, X. Jiang, F. Jia, N. Li, A. K. Nandi, Applications of machine learning to machine fault diagnosis: A review and roadmap, Mechanical Systems and Signal Processing 138 (2020) 106587.
  * [5] J. Cervantes, F. Garcia-Lamont, L. Rodríguez-Mazahua, A. Lopez, A comprehensive survey on support vector machine classification: Applications, challenges and trends, Neurocomputing 408 (2020) 189–215.
  * [6] G. A. Susto, A. Schirru, S. Pampuri, S. McLoone, A. Beghi, Machine learning for predictive maintenance: A multiple classifier approach, IEEE Transactions on Industrial Informatics 11 (3) (2014) 812–820.
  * [7] W. Wang, A model to predict the residual life of rolling element bearings given monitored condition information to date, IMA Journal of management mathematics 13 (1) (2002) 3–16.
  * [8] X. Jin, Y. Sun, Z. Que, Y. Wang, T. W. Chow, Anomaly detection and fault prognosis for bearings, IEEE Transactions on Instrumentation and Measurement 65 (9) (2016) 2046–2054.
  * [9] Y. Wang, Y. Peng, Y. Zi, X. Jin, K.-L. Tsui, A two-stage data-driven-based prognostic approach for bearing degradation problem, IEEE Transactions on industrial informatics 12 (3) (2016) 924–932.
  * [10] D. Fernandez-Francos, D. Martinez-Rego, O. Fontenla-Romero, A. Alonso-Betanzos, Automatic bearing fault diagnosis based on one-class $\nu$-svm, Computers & Industrial Engineering 64 (1) (2013) 357–365.
  * [11] S. Mahadevan, S. L. Shah, Fault detection and diagnosis in process data using one-class support vector machines, Journal of process control 19 (10) (2009) 1627–1639.
  * [12] K. K. Reddy, S. Sarkar, V. Venugopalan, M. Giering, Anomaly detection and fault disambiguation in large flight data: a multi-modal deep auto-encoder approach, in: Annual Conference of the PHM Society, Vol. 8, 2016.
  * [13] A. Heng, S. Zhang, A. C. Tan, J. Mathew, Rotating machinery prognostics: State of the art, challenges and opportunities, Mechanical systems and signal processing 23 (3) (2009) 724–739.
  * [14] J. Sikorska, M. Hodkiewicz, L. Ma, Prognostic modelling options for remaining useful life estimation by industry, Mechanical systems and signal processing 25 (5) (2011) 1803–1836.
  * [15] J. Lee, F. Wu, W. Zhao, M. Ghaffari, L. Liao, D. Siegel, Prognostics and health management design for rotary machinery systems—reviews, methodology and applications, Mechanical systems and signal processing 42 (1-2) (2014) 314–334.
  * [16] X.-S. Si, W. Wang, C.-H. Hu, D.-H. Zhou, Remaining useful life estimation–a review on the statistical data driven approaches, European journal of operational research 213 (1) (2011) 1–14.
  * [17] Y. Lei, N. Li, L. Guo, N. Li, T. Yan, J. Lin, Machinery health prognostics: A systematic review from data acquisition to rul prediction, Mechanical Systems and Signal Processing 104 (2018) 799–834.
  * [18] A. Soualhi, K. Medjaher, N. Zerhouni, Bearing health monitoring based on hilbert–huang transform, support vector machine, and regression, IEEE Transactions on Instrumentation and Measurement 64 (1) (2014) 52–62.
  * [19] R. Khelif, B. Chebel-Morello, S. Malinowski, E. Laajili, F. Fnaiech, N. Zerhouni, Direct remaining useful life estimation based on support vector regression, IEEE Transactions on industrial electronics 64 (3) (2016) 2276–2285.
  * [20] A. Saxena, K. Goebel, D. Simon, N. Eklund, Damage propagation modeling for aircraft engine run-to-failure simulation, in: 2008 international conference on prognostics and health management, IEEE, 2008, pp. 1–9.
  * [21] X. Li, Q. Ding, J.-Q. Sun, Remaining useful life estimation in prognostics using deep convolution neural networks, Reliability Engineering & System Safety 172 (2018) 1–11.
  * [22] Y. Wu, M. Yuan, S. Dong, L. Lin, Y. Liu, Remaining useful life estimation of engineered systems using vanilla lstm neural networks, Neurocomputing 275 (2018) 167–179.
  * [23] D. Wang, K.-L. Tsui, Q. Miao, Prognostics and health management: A review of vibration based bearing and gear health indicators, IEEE Access 6 (2017) 665–676.
  * [24] V. Hamaide, F. Glineur, Unsupervised minimum redundancy maximum relevance feature selection for predictive maintenance: Application to a rotating machine, International Journal of Prognostics and Health Management 12 (2) (2021).
  * [25] B. Schölkopf, J. C. Platt, J. Shawe-Taylor, A. J. Smola, R. C. Williamson, Estimating the support of a high-dimensional distribution, Neural computation 13 (7) (2001) 1443–1471.
  * [26] C.-C. Chang, C.-J. Lin, Libsvm: A library for support vector machines, ACM transactions on intelligent systems and technology (TIST) 2 (3) (2011) 1–27.
  * [27] B. E. Boser, I. M. Guyon, V. N. Vapnik, A training algorithm for optimal margin classifiers, in: Proceedings of the fifth annual workshop on Computational learning theory, ACM, 1992, pp. 144–152.
  * [28] H. Drucker, C. J. Burges, L. Kaufman, A. J. Smola, V. Vapnik, Support vector regression machines, in: Advances in neural information processing systems, 1997, pp. 155–161.
  * [29] A. J. Smola, B. Schölkopf, A tutorial on support vector regression, Statistics and computing 14 (3) (2004) 199–222.
  * [30] W. Kleeven, M. Abs, E. Forton, S. Henrotin, Y. Jongen, V. Nuttens, Y. Paradis, E. Pearson, S. Quets, J. Van de Walle, et al., The IBA superconducting synchrocyclotron project S2C2, in: Proc. Cyclotrons, 2013, pp. 115–119.
  * [31] MathWorks, Signal features - matlab & simulink (2021).
URL https://mathworks.com/help/predmaint/ug/signal-features.html

  * [32] R. R. Schoen, T. G. Habetler, F. Kamran, R. Bartfield, Motor bearing damage detection using stator current monitoring, IEEE transactions on industry applications 31 (6) (1995) 1274–1279.

Aggregation | FP | TP | B score | F score | Final score
---|---|---|---|---|---
MA 12h | 2/39 | 10/11 | 0.54 | 0.940 | 0.906
EX 12h | 2/39 | 10/11 | 0.546 | 0.940 | 0.906
MA 24 | 2/39 | 9/11 | 0.531 | 0.919 | 0.886
EX 24h | 2/39 | 10/11 | 0.54 | 0.940 | 0.906
MA 48h | 3/39 | 11/11 | 0.464 | 0.937 | 0.84
EX 48h | 4/39 | 9/11 | 0.532 | 0.880 | 0.85
MA 5D | 2/39 | 7/11 | 0.339 | 0.863 | 0.819
EX 5D | 0/39 | 6/11 | 0.265 | 0.857 | 0.806
Table 5: Binary classification results with moving average (MA) aggregation and exponential smoothing (EX) Aggregation | FP | TP | B score | F score | Final score
---|---|---|---|---|---
MA 12h | 5/39 | 9/11 | 0.438 | 0.860 | 0.824
EX 12h | 6/39 | 8/11 | 0.441 | 0.819 | 0.787
MA 24 | 6/39 | 7/11 | 0.423 | 0.793 | 0.762
EX 24h | 6/39 | 8/11 | 0.435 | 0.819 | 0.786
MA 48h | 6/39 | 8/11 | 0.458 | 0.819 | 0.788
EX 48h | 5/39 | 8/11 | 0.48 | 0.838 | 0.808
MA 5D | 2/39 | 7/11 | 0.42 | 0.863 | 0.825
EX 5D | 0/39 | 7/11 | 0.348 | 0.897 | 0.85
Table 6: Regression ReLu results Aggregation | FP | TP | B score | F score | Final score
---|---|---|---|---|---
MA 12h | 5/39 | 8/11 | 0.47 | 0.838 | 0.807
EX 12h | 5/39 | 9/11 | 0.472 | 0.860 | 0.827
MA 24 | 5/39 | 9/11 | 0.518 | 0.860 | 0.831
EX 24h | 5/39 | 10/11 | 0.525 | 0.879 | 0.848
MA 48h | 5/39 | 8/11 | 0.455 | 0.838 | 0.805
EX 48h | 5/39 | 9/11 | 0.478 | 0.860 | 0.827
MA 5D | 5/39 | 8/11 | 0.525 | 0.838 | 0.811
EX 5D | 3/39 | 8/11 | 0.433 | 0.875 | 0.838
Table 7: Multi-class classification results Aggregation | FP | TP | B score | F score | Final score
---|---|---|---|---|---
MA 12h | 2/39 | 9/11 | 0.492 | 0.919 | 0.883
EX 12h | 2/39 | 9/11 | 0.492 | 0.919 | 0.883
MA 24 | 3/39 | 9/11 | 0.496 | 0.9 | 0.865
EX 24h | 3/39 | 9/11 | 0.489 | 0.9 | 0.864
MA 48h | 2/39 | 8/11 | 0.475 | 0.894 | 0.858
EX 48h | 4/39 | 8/11 | 0.465 | 0.857 | 0.824
MA 5D | 2/39 | 7/11 | 0.385 | 0.863 | 0.822
EX 5D | 0/39 | 8/11 | 0.297 | 0.930 | 0.875
Table 8: Univariate model results