# Proton Computed Tomography Based on Richardson – Lucy Algorithm

Ákos Sudár1,2 and Gergely Gábor Barnaföldi1
for the Bergen pCT collaboration  1Wigner Research Centre for Physics,
Institute for Particle and Nuclear Physics, Budapest, Hungary 2Budapest
University of Technology and Economics , Institute of Nuclear Techniques,
Budapest, Hungary

###### Abstract.

Objective: Proton therapy is a emerging method against cancer. One of the main
development is to increase the accuracy of the Bragg-peak position
calculation, which requires more precise relative stopping power (RSP)
measurements. An excellent choice is the application of proton computed
tomography (pCT) systems which take the images under similar conditions to
treatment as they use the same irradiation device and hadron beam for imaging
and treatment. A key aim is to develop an accurate image reconstruction
algorithm for pCT systems to reach their maximal performance.

Approach: An image reconstruction algorithm was developed in this work, which
is suitable to reconstruct pCT images from the energy, position and direction
measurement of individual protons. The flexibility of an iterative image
reconstruction algorithm was utilised to appropriately model the trajectory of
protons. Monte Carlo (MC) simulations of a Derenzo and a CTP404 phantom was
used to test the accuracy of the image reconstruction.

Main results: The Richardson – Lucy algorithm was applied first and
successfully for pCT image reconstruction. Probability density based approach
was applied for interaction (system) matrix generation, which is an advanced
approach to consider the uncertain path of the protons in the patient.

Significance: The track of protons are scattered when they travel through
material at the hadron therapy energies. This property limits the achievable
spatial resolution, especially for single-sided pCT setups investigated in
this study. The main motivation of the presented research is to test new
approaches for the image reconstruction, focusing on the achieved spatial- and
density resolution and the image noise. Realistic imaging setup were simulated
with reasonably low proton statistics, to achieve results, which is likely to
be reproducible in clinical environment.

##  1\. Introduction

Hadron therapy is an emerging and efficient curing method against cancer. The
increasing number of hadron therapy centers and the number of successful
treatments presents its success. Today’s accelerator techniques led us to use
protons or heavier ions as bombarding particles. The application of hadron
beams instead of X-ray results more focused dose distribution [6]. Indeed,
using higher mass number beams than proton (He, C and O) can results an
increased relative biological effectiveness (RBE) in the tumor volume, while
keep the RBE close to one in healthy tissues [7]. The higher the dose gradient
around the treated volume requires the lower uncertainty in the relative
stopping power (RSP) distribution during dose planning, to avoid insufficient
dosage of the tumor or the overdose of organs at risk.

The development of proton Computed Tomography (pCT) techniques are promising
solutions for the above problems. Applying the same irradiation device, beam,
and hadron for both the medical imaging and the treatments can significantly
reduce the uncertainties of the imaging. To obtain this, two main imaging
strategies are exist:

  * (i)

The first concept is to measure the average energy loss of the proton beam.
This design is feasible from a technical point of view, but can only achieve
poor spatial resolution with the clinically available proton beams [15].

  * (ii)

The second concept is the so called list mode imaging concept, which measures
the energy loss and in parallel it estimates the path of each individual
proton. Monte Carlo (MC) simulations and prototype measurements showed that
this solution can meet the required spatial and density resolutions, so the
focus is moved toward this direction [13].

Nowadays, the pCT scanner R&Ds around the world are tending to reach the
prototyping and clinical/pre-clinical testing phase, which requires the
integration of the prototype scanners into clinical environment. Following the
list mode strategy, the path estimation of individual protons is usually based
on the measurements of the upstream and downstream tracker detector pairs,
which concept is called the double-sided scanner design (figure 1). One
important further step can be the abandonment of the upstream tracker
detectors and the application of a single-sided scanner design. The drawback
of this latter concept is the less accurate proton path measurement, however
the study by [27] concluded that the achievable spatial resolution meets the
minimum requirement. Nevertheless with lower-precision proton path measurement
compromises the spatial resolution of the scanner which immediately motivates
the development and application of more accurate image reconstruction
algorithms.

Figure 1. Detector design of the list mode imaging concept.

The realistic clinical applicability requires to complete the data taking in
minutes, which result in 1-10 million protons per second measurement rate. The
ultimate goal would be to finish the image capturing within the minimum gantry
rotation time. The LLU/UCSC Phase-II Scanner prototype detector showed speed
up to 1.2 million proton per second, which probably can be upgraded with 50%
in the near future [14]. To increase the data taking rate to 10 million
protons per second, two possible directions exist: the first is to apply
faster readout frequency of 10 MHz at least, the second is to measure multiple
proton tracks within one readout frame. The second solution fits mostly with
the single-sided scanner design, because this solution avoids the pairing
problem of the upstream and downstream measurements, which leads to track
confusion in case of a double-sided scanner even with low number of protons in
a frame.

Multiple proton measurements fits best for silicon pixel trackers and silicon
pixel sensor based range counters as presented by the Bergen pCT Collaboration
[19, 20, 3]. However, multiple proton measurements can be done by applying
three silicon strip detectors rotated relative to each other as presented by
the PRaVDA collaboration [8]. Another layout has been designed by the iMPACT
group. The ProXY detector combines the two acceleration possibilities with
monolithic active pixel detectors. Applying this layout, about 50 MHz readout
frequency can be reached, and it is planned to measure multiple-proton events
[17].

In this paper authors present for the first time a novel points, which may be
applied in the proposed pCT detector concepts: the Richardson – Lucy algorithm
[22, 16] was applied for the image reconstruction of a single-sided pCT
scanner. The paper organized as follows: section 2 begins with the general
approach to the image reconstruction problem itself, followed by the
presentation of the details of the Richardson – Lucy algorithm and the proton-
phantom interaction model. Section 3 compares detector designs and presents
the applied Monte Carlo simulation. Section 3 also contains the evaluation of
the spatial- and density resolution of phantoms. Results are summarized and
discussed in section 4 and 5, respectively.

##  2\. The Image Reconstruction Algorithm

The role of the image reconstruction is to give back the relative stopping
power (RSP) distribution from the measured data. Two family of image
reconstruction techniques exist: the first family contains the filtered
backprojections, in contrast with the second includes the iterative
reconstructions. The first family usually use integrals along straight lines,
which seems to be an inaccurate approximation for the scattered proton
trajectory. The so-called distance-driven backprojection belongs to this
family but can take into account the curvature of the MLPs during the filtered
backprojection [23]. This method provides reasonably good spatial resolution
however, it requires very high statistics, which might be not accepted by
requirements of the clinical use. The second family of the image
reconstructions models the imaging as the interaction of proton tracks and
volumetric pixels (voxels) of the reconstruction space. This approach is
suitable to handle curved proton trajectories and reaches reasonably good
spatial and density resolution with acceptable statistics, however requires
higher computational power. This method models the imaging as a large linear
system of equations, which is described by the following general algebraic
form:

$\textbf{y}=\textbf{A}\cdot\textbf{x},$ (1)

where y is an $m$-dimensional vector, which has typically $10^{8}$-$10^{9}$
elements. The y contains the water equivalent path length (WEPL) reduction of
the protons in the reconstruction area. Variable x is an $n$-dimensional
vector (typically $10^{5}$-$10^{7}$ elements) containing the relative stopping
power (RSP) of the voxels. Finally A is the so called system matrix, which has
$n\times m$ elements of $10^{13}$-$10^{16}$. The system matrix contains the
interaction coefficients between protons and voxels. The information in matrix
A can be described as the (expected) length of the proton’s path in the voxel.
In practice, $m$ is usually larger than $n$, so the linear equation system is
over-determined. The goal of the image reconstruction in general, is to
determinte the values of vector, x with the knowledge of vector, y and system
matrix, A.

Orthogonal projection based iterative algorithms are widely used for pCT image
reconstruction as presented by [13, 9, 18, 4, 10]. In this work the authors
applied the Richardson – Lucy deconvolution, which reached reasonable quality
increasement in the field of emission tomography [26], and have not been used
earlier for proton computed tomography.

###  2.1. The Richardson – Lucy Algorithm

The Richardson – Lucy deconvolution iteration cycle [22, 16] originating from
the field of optics, and known as a fixed point iteration. The iterative
solution is based on the formula,

$\textbf{x}_{i}^{k+1}=\textbf{x}_{i}^{k}\frac{1}{\sum\limits_{j}\textbf{A}_{ij}}\sum\limits_{j}\frac{\textbf{y}_{j}}{\sum\limits_{l}\textbf{A}_{lj}\,\textbf{x}_{l}^{k}}\textbf{A}_{ij}~{},$
(2)

for every $i=1,\ ...,\ N$, where $N$ is the length of vector x, which contains
the RSP of the voxels, $k$ is the number of iterations, matrix
$\textbf{A}_{ij}$ contains the interaction coefficients between the proton
trajectories and the voxels, $j=1,\ ...,\ M$ is the index of the trajectories,
where $M$ is the number of the trajectories, $\textbf{y}_{j}$ contains the
integrated RSP along the trajectories, which is equivalent with the WEPL
reduction of the protons travelling along the trajectories. The
$\textbf{y}_{j}/\sum\limits_{l}\textbf{A}_{lj}\,\textbf{x}_{l}^{k}$ term is
usually called Hadamard ratio, and it represents the ratio of the integrated
RSP along the proton path and its estimate based on the voxel values
calculated in the previous iteration.

###  2.2. The Proton-Phantom Interaction

Instead of the most simple straight line approximation, the literature uses
the estimated (most likely) path of the protons, based on the upstream and
downstream measurements of proton track position and angle in case of a
double-sided scanner design. Novel formulae are available in [24, 30, 25, 15]
to calculate the most likely path (MLP) of the protons. In case of a single-
sided scanner design, where upstream measurements are not available, the beam
information is used. Certainly, this contains much more uncertainty than a
precise measurement. In this work the authors followed the formalism of [15]
as it considers the uncertainty of the measurements and the beam. The path of
the protons was considered to be straight outside the phantom.

In this article the authors applied an advanced approach suggested by [30]: a
Gaussian probability density distribution of the real proton path around the
MLP, which takes into account the uncertain path of the protons in the
phantom. This approach was used by [29] for pCT image reconstruction. An
average standard deviation was considered during the proton path in the
patient, which is an approximation compared to the depth dependent probability
density investigated by [30], which work did not deal with image
reconstruction. The average standard deviation was chosen based on the
experience of the authors during the development.

##  3\. Simulations with the Algorithm

The ultimate goal of pCT imaging for proton therapy is to provide a stable
basis for accurate dose planning. However, it is a challenge to define a
measure, which characterizes the goodness of a reconstructed RSP distribution
for dose calculation. Instead of this missing ideal measurement, general image
properties (spatial and density resolution, image noise) are usually used to
quantify image quality. To study this, the imaging of dedicated spatial and
density resolution phantoms was simulated with Monte Carlo techniques,
reconstructed by the formerly described method and evaluated following the
instructions later on this section.

###  3.1. The Proton CT Scanner Model

A single-sided detector design (figure 2) with a 230 MeV/u pencil beam was
investigated. The full width at half maximum (FWHM) of the Gaussian beam was 7
mm (with about 3 mm standard deviation), the spot divergence was set to 2.8
mrad and the spot emittance was 3.0 mrad$\times$mm, following the beam model
of [27]. Three different detector layer was investigated: the first is an
idealized detector with no measurement errors, the second is a silicon pixel
tracker modelled after the design of the Bergen pCT Collaboration [20, 3] and
the third is a silicon strip detector based tracker layer followed the
LLU/UCSC Phase-II Scanner design of the Loma Linda University (LLU) and the
University of California at Santa Cruz (UCSC) [14]. We note, that the results
of this work is valid for an envisioned single-sided scanner, built of the
LLU/UCSC Phase-II Scanner tracker layers, in comparison with the existing
LLU/UCSC Phase-II Scanner which is a state-of-the-art double-sided setup.

Figure 2. Single-sided list mode detector design.

The properties of the three detector layer setups, the idealized, the silicon
pixel, and the silicon strip, are summarized in table 1. The idealized setup
is a single sensitive layer, but in the latter two realistic cases each
tracker layer contains two sensitive planes due to the existing technological
solutions (figure 2). If the detection is based on a silicon pixel detector a
double structure of two equivalent sensitive planes need to apply to fully
cover the alternating sensitive and readout electronics panels. While applying
silicon strip detectors, two separate planes are required for the perpenicular
$x$ and $y$ directions. Table 1 contains the joint material budget of these
double layers. The WEPL resolution of both realistic setups was chosen to be 3
mm (standard deviation of a normal distribution), which was added to the
simulated range straggling in the phantom. This is a realistic uncertainty,
however this is a very rudimentary model, as the measurement error is likely
to depend on the remaining range of the protons behind the patient. The
distance of the first detector pair to the rotation axis (isocenter) was
chosen for 400 mm in all cases, which result in 300 mm and 325 mm detector-
phantom distance for the Derenzo and the CTP404 phantoms, respectively.
Similar distances are used for portal detectors applied in photon therapy
gantries, so they can be considered to be realistic for the future pCT devices
as well [15].

| Unit | Ideal setup | Silicon pixel | Silicon strip
---|---|---|---|---
Layer material budget ($x/X_{0}$) | - | 0 | $4.2\times 10^{-3}$ | $8.5\times 10^{-3}$
Distance between layers | mm | - | 50 | 50
Spatial resolution |  $\mu$m | 0 | 5 | 66
Angular resolution (130-230 MeV/u) | mrad | 0 | 1.7-2.9 | 3.1-4.6
Correlation (130-230 MeV/u) |  mrad$\times$mm | 0 | $-5\times 10^{-4}$ | $-8.7\times 10^{-2}$
Statistical WEPL resolution | mm | 0 | 3.0 | 3.0
Table 1. Comparison of tracker detector pair model parameters: Ideal setup,
with no measurements errors, Silicon pixel detector based on the design of the
Bergen pCT Collaboration [20, 3], and Silicon strip detector model following
the structure of the LLU/UCSC Phase-II Scanner [14].

###  3.2. The Applied Phantoms

To test and validate the application of Richardson – Lucy algorithm was
required to apply standardised evaluation methods. Therefore, we applied two
widely-applied phantom in our analysis. In this study the RSP distribution was
reconstructed in one plane of the phantoms, so the phantoms were considered to
be offset invariant in the direction of the rotation axis. To ensure the
offset invariance 400 mm high phantoms were simulated in the axis direction.

The spatial resolution of the reconstruction was measured with the MC imaging
of the Derenzo phantom: a 200 mm diameter water cylinder, which contains six
sectors of 1.5-6 mm diameters aluminium rods, specially chosen for the current
analysis. The original idea of this phantom comes by [5].

We also used a CTP404 phantom for our study. CTP404 is produced by [28], and
designed to measure how accurately a material property is reconstructed in a
homogeneous region of the phantom. The reconstruction accuracy of the RSP can
be evaluated for proton CT imaging, also referred to as density resolution in
the literature. The CTP404 phantom is an epoxy 150 mm diameter cylinder, which
contains 8 different material inserts with a diameter of 12.2 mm. The average
RSP of the inserts was evaluated in an 8 mm diameter circle in the middle of
the inserts and compared to the real RSP values investigated by [3]. The
standard deviation of the RSP was also evaluated in every inserts to
characterise the noise of the image.

###  3.3. Steps of the Simulation with the Algorithm

A simulation code was developed to test the Richardson – Lucy algorithm, which
was divided into the following steps (schematic in figure 3):

Figure 3. Simulation steps.

  * (1)

The data taking was simulated with Monte Carlo method. The beam and the
phantom were modeled appropriately in the simulation. The Geant4 (version
11.0.0) [1, 2] was used with GATE (version 9.2) [12, 11]. In the reference
physics list settings QGSP_BIC_EMY was activated for the calculations. Data
taking of one slice was simulated from 180 directions in $2^{\circ}$ steps.
Field of view (FOV) with 220 mm was applied for the Derenzo phantom (111 beam
positions with 2 mm steps), and 170 mm FOV was used in case of the CTP404
phantom (86 beam positions with 2 mm steps). Overal $\sim 2$ million and $\sim
1.5$ million primary protons were simulated from which $\sim 1.2$ million and
$\sim 0.95$ million remains after 3 sigma filtering in case of the Derenzo and
the CTP404 phantoms, respectively. Instead of modeling the detector in the MC
simulation, the exact position, direction and energy of the protons were read
out at the position of the first tracker layer, and the measurement
uncertainties were assigned in the next step.

  * (2)

In this step the errors of position and direction measurements were drawn from
correlated Gaussian distributions and added to the exact positions and
directions of simulated protons. The measurement uncertainty was calculated
based on the guideline of [15]. The WEPL measurement error also was randomly
assigned from a Gaussian distribution to the WEPL of the protons calculated
from their energy losses, simulated in the previous step. The parameters from
table 1 were used. In case of the ideal setup this step was certainly skipped,
since the lack of errors/uncertainties of the idealized case.

  * (3)

A 3 sigma filtering is applied for the direction and WEPL of the protons
originating from the same beam spot. The goal of this step was to filter out
protons, which undergo nuclear collisions in the patient. This type of
filtering was suggested and used by [25].

  * (4)

Calculation of the most probable incoming and outgoing position of the protons
on a cylinder around the phantom is performed. In this step the formalism of
[15] was applied. The diameter of the cylinder was chosen to be 10 mm wider
than that of the phantom in order to avoid artefacts.

  * (5)

The Richardson – Lucy algorithm was used to reconstruct the RSP distribution
from the individual proton histories. On the fly system matrix calculation was
applied based on simplified probability density around a third order spline
approximation of the MLP.

  * (6)

In the final step the spatial resolution was evaluated based on the
reconstruction of the Derenzo phantom. The density resolution and the image
noise were calculated from the reconstructed CTP404 phantom.

Calculations were done on the machines of the Wigner Scientific Computing
Laboratory’s hardware. The computationally demanding part of the algorithm was
running on four 1080 Ti GPU cards. The focus of the current work was on the
proof of the concept of probability density based system matrix calculation
with Richardson – Lucy reconstruction algorithm, so the authors did not focus
on the optimization of the implemented code.

###  3.4. Evalution of the Derenzo Phantom

The evaluation of the measured phantom is based on the above simulation, and
starts with the comparison of the reconstructed intensity in the position of
the aluminium rods (peaks) and medium between them (valleys) as it is
demonstrated in figure 4. After subtraction of the background (defined by the
tails of the distribution) the valley-to-peak ratio can be calculated.

Figure 4. The RSP distribution along the sidelines (with three different
colors) of the triangle from the 4 mm rods in the reconstructed Derenzo
phantom. The Bergen pCT setup was applied.

The blurring effect of the reconstruction is modeled as a convolution with the
so called point-spread function. It is a Gaussian function, and its Fourier
transform is the modulation transfer function (MTF). The frequency of the MTF
at 10 % (measured in units of line pair per cm) can be derived from the
valley-to-peak ratio (figure 5) and quantifies the spatial resolution of the
reconstructed image. If the valley-to-peak ratio is too close to zero or one
the image noise suppress the information about the point-spread function, so a
sector of the Derenzo phantom is suitable to cover only a limited resolution
range. The phantom contains six sectors with different rod diameters to
increase the range of resolution that can be evaluated.

Figure 5. The spatial resolution (in the units of linepair per cm) as a
function of the valley to peak ratio.

##  4\. Results

As we pointed our earlier, in our study we investigated the simplified model
with a single-sided detector setup. The center line of all beams falls into
one perpendicular plane to the rotation axis. This image slice was
reconstructed containing 256$\times$256 pixels with 1 mm2 size. Every proton
was assigned to this layer, without take into account the deviation of their
path in the direction of the rotation axis. In the reconstruction step, every
reconstructed image was evaluated after 600 iterations as an optimum between
the magnitude of spatial resolution and the image noise. The standard
deviation ($\sigma$) around the MLP of the protons were set about 0.4 mm, 0.5
mm and 1.0 mm according to the spatial uncertainty of the proton path for the
ideal, silicon pixel and silicon strip detector layers, respectively.

The result of the reconstructed images of the Derenzo and CTP404 phantoms are
shown in figure 6, on the left and right columns, respectively.

Figure 6. Richardson – Lucy algorithm based reconstruction of the Derenzo and
CTP404 phantoms is presented on the left and right columns, respectively. From
top to bottom the ideal, the pixel detector and the strip detector layers were
drawn.

Top row of the figure 6 presentes the idealized case, which certainly the most
clearest reconstruction for both Derenzo and CTP404 phantoms. The middle and
bottom rows are present more realistic imaging models with silicon pixel and
silicon strip detectors, respectively. It is clearly visible that as the
position and direction measurement uncertainties increases, the spatial
resolution becomes worse and worse.

To quantify the quality of the above reconstructed images, the iteration-by-
iteration evolution of the spatial resolution, the image noise and the density
resolution are drawn in figure 7, respectively, form left to right. The
spatial resolution (left panel) of the ideal setup was found to be 2.4 lp/cm
for the ideal setup, 2.0 lp/cm and 1.5 lp/cm for the silicon pixel and the
silicon strip detector based setups, respectively. The noise (middle panel) of
the ideal and silicon pixel setups seems to be similar (5 % and 4.8 %,
respectively), while the silicon strip detector based setup reached
significantly lower noise (3.3 %). This observation indicates that if the
number of protons are fixed, then the $\sigma$ of the reconstruction has the
most important effect on noise, instead of spatial or WEPL measurement
uncertainties. We note, that the value of $\sigma$ parameter was set
significantly higher in case of the silicon strip setup, following the higher
uncertainty in the proton path measurement compared with the other two
detector layers.

Figure 7. Left: the spatial resolution, middle: the image noise and right: the
average relative RSP error as a function of the iteration number.

The error of the density reconstruction is quickly decreasing up to 70-90
iterations, thereafter all saturates, while for the ideal setup it is getting
better at large iteration numbers. The average relative RSP difference (except
air) was found to be 0.3 % for the ideal and 0.5 % for the realistic setups
after 600 iterations. The RSP of the air instead of the real 0.001 was found
to be 0.036, 0.051 and 0.061 for the ideal, the silicon pixel and the silicon
strip setups. The density resolution of the tissues around and above water can
be obviously reconstructed with better accuracy than the required 1 % [21],
but the density resolution is pure for low density regions. The reconstructed
RSP of the air inserts was significantly decreasing even after 600 iterations,
so even more iterations may solve this uncertainty. The RSP of all inserts
(instead of air) was underestimated, most likely caused by the overestimated
RSP of the air around the phantom. If the reconstruction area would be limited
to the area of the phantom (instead of the whole $256\times 256$ mm2), the RSP
estimation of these inserts probably would be even more accurate.

##  5\. Discussion

Based on our simulation, we found that the Richardson – Lucy algorithm with
probability density based proton-phantom interaction reached 2.4 lp/cm spatial
resolution for the ideal and 2.0 lp/cm resolution for a realistic setup. The
density resolution was found to be significantly better than the required 1 %
RSP accuracy, except of the low density region. We observed, that, the number
of iterations further improves the spatial resolutions and in parallel the
density resolution at low RSP as well. Meanwhile, the larger noise (around 5 %
after 600 iterations) limits the applicable number of iterations. This
limitation could be exceeded by imaging with higher statistics or maybe by
implementing superiorization, which reduces the noise.

##  6\. Summary

In this work the application of the Richardson – Lucy algorithm with
probability density based proton-phantom interaction calculation for proton CT
image reconstruction has been presented. The authors applied clinically
realistic setups and parameters in the Monte Carlo simulations: 400 mm
detector-isocenter distance, 1.5-2 million primary protons per image slice and
realistic beam and detector characteristics. For testing and for the
evaluation of the resolutions two widely used phantom were applied the Derenzo
and the CTP404.

The authors concluded that the presented reconstruction method meets the
required density resolution, and have similar density resolution as the state
of art prototypes. However, the spatial resolution of the images are
promising, but has not reached the clinical requirements yet. Further
development of the algorithm is necessary. It is also important to mention,
that however the reconstructed images are noisy due to the low statistics,
they are almost artefact free without the application of any post processing
method.

The authors are continuing the algorithm development, with a focus on the
spatial resolution and reconstruction time, which limited the investigations
of the current work to only one layer. Indeed the possibility of the speedup
of the algorithm has been also planned.

## Members of the Bergen pCT Collaboration

Max Aehlea, Johan Almeb, Gergely Gábor Barnaföldic, Tea Bodovab, Vyacheslav
Borshchovd, Anthony van den Brinke, Mamdouh Chaarb, Viljar Eikelande, Gregory
Feofilovf, Christoph Garthg, Nicolas R. Gaugera, Georgi Genovb, Ola Grøttvikb,
Havard Helstruph, Sergey Igolkinf, Ralf Keideli, Chinorat Kobdajj, Tobias
Kortusi, Viktor Leonhardtg, Shruti Mehendaleb, Raju Ningappa Mulawadei, Odd
Harald Odlandk, b, George O’Neillb, Gábor Pappl, Thomas Peitzmanne, Helge Egil
Seime Pettersenk, Pierluigi Piersimonib,m, Maksym Protsenkod, Max Rauchb,
Attiq Ur Rehmanb, Matthias Richtern, Dieter Röhrichb, Joshua Santanai,
Alexander Schillingi, Joao Secoo, p, Arnon Songmoolnakb, j, Jarle Rambo
Sølieq, Ákos Sudárc, r Ganesh Tambaveb, Ihor Tymchukd, Kjetil Ullalandb,
Monika Varga-Kofaragoc, Lennart Volzs, t, Boris Wagnerb, Steffen Wendzeli,
Alexander Wiebeli, RenZheng Xiaob, u, Shiming Yangb, Hiroki Yokoyamae,
Sebastian Zillieni

a) Chair for Scientific Computing, TU Kaiserslautern, 67663 Kaiserslautern,
Germany; b) Department of Physics and Technology, University of Bergen, 5007
Bergen, Norway; c) Wigner Research Centre for Physics, Budapest, Hungary; d)
Research and Production Enterprise “LTU” (RPELTU), Kharkiv, Ukraine; e)
Institute for Subatomic Physics, Utrecht University/Nikhef, Utrecht,
Netherlands; f) St. Petersburg University, St. Petersburg, Russia; g)
Scientific Visualization Lab, TU Kaiserslautern, 67663 Kaiserslautern,
Germany; h) Department of Computer Science, Electrical Engineering and
Mathematical Sciences, Western Norway University of Applied Sciences, 5020
Bergen, Norway; i) Center for Technology and Transfer (ZTT), University of
Applied Sciences Worms, Worms, Germany; j) Institute of Science, Suranaree
University of Technology, Nakhon Ratchasima, Thailand; k) Department of
Oncology and Medical Physics, Haukeland University Hospital, 5021 Bergen,
Norway; l) Institute for Physics, Eötvös Lóránd University, 1/A Pázmány P.
Sétány, H-1117 Budapest, Hungary; m) UniCamillus – Saint Camillus
International University of Health Sciences, Rome, Italy; n) Department of
Physics, University of Oslo, 0371 Oslo, Norway; o) Department of Biomedical
Physics in Radiation Oncology, DKFZ—German Cancer Research Center, Heidelberg,
Germany; p) Department of Physics and Astronomy, Heidelberg University,
Heidelberg, Germany; q) Department of Diagnostic Physics, Division of
Radiology and Nuclear Medicine, Oslo University Hospital, Oslo, Norway; r)
Budapest University of Technology and Economics, Budapest, Hungary; s)
Biophysics, GSI Helmholtz Center for Heavy Ion Research GmbH, Darmstadt,
Germany; t) Department of Medical Physics and Biomedical Engineering,
University College London, London, UK; u) College of Mechanical & Power
Engineering, China Three Gorges University, Yichang, People’s Republic of
China

## Acknowledgement

The authors would like to thank the support of the Hungarian National
Research, Development and Innovation Office (NKFIH) grants under the contract
numbers OTKA K135515 and 2019-2.1.6-NEMZ_KI-2019-00011,
2020-2.1.1-ED-2021-00179. This work was also supported by the Research Council
of Norway (Norges forskningsrad) and the University of Bergen, grant number
250858. The authors acknowledge the support os Trond Mohn Foundation
(BFS2017TMT07). Computational resources were provided by the Wigner Scientific
Computing Laboratory (WSCLAB).

## References

  * [1] S. Agostinelli et al.  “Geant4—a simulation toolkit”  In _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 506.3, 2003, pp. 250–303  DOI: https://doi.org/10.1016/S0168-9002(03)01368-8
  * [2] J. Allison et al.  “Geant4 developments and applications”  In _IEEE Transactions on Nuclear Science_ 53.1, 2006, pp. 270–278  DOI: 10.1109/TNS.2006.869826
  * [3] Johan Alme et al.  “A High-Granularity Digital Tracking Calorimeter Optimized for Proton CT”  In _Frontiers in Physics_ 8, 2020  DOI: 10.3389/fphy.2020.568243
  * [4] Yair Censor, Tommy Elfving, Gabor T. Herman and Touraj Nikazad  “On Diagonally Relaxed Orthogonal Projection Methods”  In _SIAM Journal on Scientific Computing_ 30, 2008  DOI: 10.1137/050639399
  * [5] S.. Derenzo et al.  “High Resolution Computed Tomography of Positron Emitters”  In _IEEE Transactions on Nuclear Science_ 24.1, 1977, pp. 544–558  DOI: 10.1109/TNS.1977.4328738
  * [6] M. Durante, R. Orecchia and JS. Loeffler  “Charged-particle therapy in cancer: clinical uses and future perspectives”  In _Nat Rev Clin Oncol_ 14.8, 2017, pp. 483–495  DOI: 10.1038/nrclinonc.2017.30
  * [7] Marco Durante, Jürgen Debus and Jay S. Loeffler  “Physics and biomedical challenges of cancer therapy with accelerated heavy ions”  In _Nature Reviews Physics_ 3.12, 2021, pp. 777–790  DOI: 10.1038/s42254-021-00368-5
  * [8] Michela Esposito et al.  “PRaVDA: The first solid-state system for proton computed tomography”  In _Physica Medica_ 55, 2018, pp. 149–154  DOI: https://doi.org/10.1016/j.ejmp.2018.10.020
  * [9] Richard Gordon, Robert Bender and Gabor T. Herman  “Algebraic Reconstruction Techniques (ART) for three-dimensional electron microscopy and X-ray photography”  In _Journal of Theoretical Biology_ 29.3, 1970, pp. 471–481  DOI: https://doi.org/10.1016/0022-5193(70)90109-8
  * [10] G.. Herman  “Fundamentals of Computerized Tomography”  Springer, 2009
  * [11] S. Jan et al.  “GATE V6: a major enhancement of the GATE simulation platform enabling modelling of CT and radiotherapy”  In _Physics in Medicine and Biology 2011-jan 20 vol. 56 iss. 4_ 56, 2011  DOI: 10.1088/0031-9155/56/4/001
  * [12] S. Jan et al.  “GATE: a simulation toolkit for PET and SPECT”  In _Physics in Medicine and Biology 2004-sep 10 vol. 49 iss. 19_ 49, 2004  DOI: 10.1088/0031-9155/49/19/007
  * [13] Robert P Johnson  “Review of medical radiography and tomography with proton beams”  In _Reports on Progress in Physics_ 81, 2017  DOI: 10.1088/1361-6633/aa8b1d
  * [14] Robert P. Johnson et al.  “A Fast Experimental Scanner for Proton CT: Technical Performance and First Experience With Phantom Scans”  In _IEEE Transactions on Nuclear Science_ 63, 2016  DOI: 10.1109/TNS.2015.2491918
  * [15] Nils Krah et al.  “A comprehensive theoretical comparison of proton imaging set-ups in terms of spatial resolution”  In _Physics in Medicine and Biology_ 63, 2018  DOI: 10.1088/1361-6560/aaca1f
  * [16] L.. Lucy  “An iterative technique for the rectification of observed distributions” Provided by the SAO/NASA Astrophysics Data System  In _Astronomical Journal_ 79, 1974, pp. 745  DOI: 10.1086/111605
  * [17] S. Mattiazzo et al.  “Advanced proton imaging in computed tomography”  In _Radiation Protection Dosimetry_ 166, 2015  DOI: 10.1093/rpd/ncv197
  * [18] S.. Penfold, R.. Schulte, Y. Censor and A.. Rosenfeld  “Total variation superiorization schemes in proton computed tomography image reconstruction”  In _Medical Physics_ 37, 2010  DOI: 10.1118/1.3504603
  * [19] H.E.S. Pettersen et al.  “Proton tracking in a high-granularity Digital Tracking Calorimeter for proton CT purposes”  In _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 860, 2017, pp. 51–61  DOI: https://doi.org/10.1016/j.nima.2017.02.007
  * [20] Helge Egil Seime Pettersen et al.  “Design optimization of a pixel-based range telescope for proton computed tomography”  In _Physica Medica_ 63, 2019, pp. 87–97  DOI: https://doi.org/10.1016/j.ejmp.2019.05.026
  * [21] G Poludniowski, N M Allinson and P M Evans  “Proton radiography and tomography with application to proton therapy”  In _British Journal of Radiology_ 88, 2015  DOI: 10.1259/bjr.20150134
  * [22] William Hadley Richardson  “Bayesian-Based Iterative Method of Image Restoration”  In _Journal of the Optical Society of America_ 62, 1972  DOI: 10.1364/JOSA.62.000055
  * [23] Simon Rit et al.  “Filtered backprojection proton CT reconstruction along most likely paths”  In _Medical Physics_ 40, 2013  DOI: 10.1118/1.4789589
  * [24] Uwe Schneider  “Multiple Coulomb scattering and spatial resolution in proton radiography”  In _Medical Physics_ 21, 1994  DOI: 10.1118/1.597212
  * [25] R.. Schulte, S.. Penfold, J.. Tafas and K.. Schubert  “A maximum likelihood proton path formalism for application in proton computed tomography”  In _Medical Physics_ 35, 2008  DOI: 10.1118/1.2986139
  * [26] L.. Shepp and Y. Vardi  “Maximum Likelihood Reconstruction for Emission Tomography”  In _IEEE Transactions on Medical Imaging_ 1.2, 1982, pp. 113–122  DOI: 10.1109/TMI.1982.4307558
  * [27] Jarle Rambo Sølie et al.  “Image quality of list-mode proton imaging without front trackers”  In _Physics in Medicine and Biology_ 65.13 IOP Publishing, 2020, pp. 135012  DOI: 10.1088/1361-6560/ab8ddb
  * [28] The Phantom Laboratory  “Catphan® 600 phantom (containing CTP404 as a section)” Last accessed 5 October 2022, 2022  URL: https://www.phantomlab.com/catphan-600
  * [29] Dongxu Wang, T. Mackie and Wolfgang A. Tomé  “On the use of a proton path probability map for proton computed tomography reconstruction”  In _Medical physics_ 37, 2010  DOI: 10.1118/1.3453767
  * [30] D C Williams  “The most likely path of an energetic charged particle through a uniform medium”  In _Physics in Medicine and Biology_ 49, 2004  DOI: 10.1088/0031-9155/49/13/010


#  Deep Compression of Neural Networks for Fault Detection on Tennessee
Eastman Chemical Processes

Mingxuan Li
Department of Computer Science
University of North Carolina at Chapel Hill
Chapel Hill, NC, USA
<EMAIL_ADDRESS>* Yuanxun Shao  Nuro Inc.
Mountain View, CA, USA
<EMAIL_ADDRESS>
###### Abstract

Artificial neural network has achieved the state-of-art performance in fault
detection on the Tennessee Eastman process, but it often requires enormous
memory to fund its massive parameters. In order to implement online real-time
fault detection, three deep compression techniques (pruning, clustering, and
quantization) are applied to reduce the computational burden. We have
extensively studied 7 different combinations of compression techniques, all
methods achieve high model compression rates over 64% while maintain high
fault detection accuracy. The best result is applying all three techniques,
which reduces the model sizes by 91.5% and remains a high accuracy over 94%.
This result leads to a smaller storage requirement in production environments,
and makes the deployment smoother in real world.

###### Index Terms:

Deep Compression, Neural Networks, Fault Detection

##  I Introduction

Artificial Neural Network (ANN) is a powerful technique that can classify
input data into predefined categories. ANN has been successfully deployed in
an enormous amount of real-life applications with superior results, including
computer vision [1], natural language processing [2], autonomous driving [3].
For a classification problem, input data is processed through a series of
weight matrices, bias vectors, and nonlinear activation functions. Following
by a classifier, ANN eventually calculates a likelihood score for each
category. During training, ANN maximizes the likelihood of the true category
(or minimize the negative of this likelihood) for taring data sets via
variants of stochastic gradient decent methods. In terms of the inference, we
run the forward path and classify testing data to the category with the
highest probability score. Even though ANN has shown great success in many
fields, sometimes the enormous model size caused by a large number of weights
can result in overwhelming computational burden (e.g. taking huge memory and
being very power demanding). In certain application scenarios that have memory
or energy consumption limitations, such as mobile platform or browser-based
systems, a smaller size ANN with high accuracy is not only necessary but also
hold in high demand.

Deep compression is a powerful technique for reducing the size of an ANN while
remaining high accuracy. There are three techniques can be used to achieve
this goal: pruning, clustering, and quantization. Pruning is the process of
removing relatively insignificant weights from weight matrices. Since the
matrices can be stored in the sparse matrix data type, the removal of any item
in one matrix directly reduces its size. Clustering, also known as weight
sharing, is the process of using the same value for multiple weights. Since
one reference can be used for several weights using this technique, it can
reduce the storage space for storing weights. At last, quantization is the
technique of storing weights in data types with lower precision in exchange
for smaller space in storage. By combing those three techniques, we are able
to reduce the size of ANN while maintaining a similar accuracy.

Fault detection, as an essential step of industrial production to ensure
efficiency and safety, has been an active field of research for past years. As
the complexity of production increases dramatically due to more advanced
technologies, deep learning has become an increasingly popular candidate for
fault detection applications due to its potential in handling complex systems,
including a large amount of measurements and many failure types for processes.
In this paper, we focus on the Tennessee Eastman process, a representative
case for fault detection and diagnostic. The Tennessee Eastman process is a
typical industrial process that consists of five main process units. It has 52
measurements and 21 fault types as defined by Downs and Vogel in [4].
Moreover, as a widely used case for fault detection, the usage of ANN to
diagnose fault has been proved to have satisfying results. To achieve the best
accuracy, instead of using a single ANN to classify the fault type, the result
in [5] uses different ANNs for different fault types and achieves 97.73%
average accuracy. However, this method requires one deep learning model for
each fault type, which requires large storage space and computational efforts
in total. Tennessee Eastman process requires the online detection ability with
low latency. Thus, the ANN pruning is an ideal candidate for reducing the
computational burden, speeding up the online inference time, as well as
maintaining the classification accuracy.

In this paper, deep compression of artificial neural networks (ANN) is used
for the fault detection on Tennessee Eastman chemical process to enable faster
online inferences with high accuracy. Sections II and III introduce the basics
of deep learning and deep compression. Section IV discusses the details of
deep learning architecture and deep compression results for fault detection of
TN Eastman chemical processes. Finally, Section V summaries this paper.

##  II deep learning

An ANN usually has 4 different types of layer: input layer, hidden layer,
softmax layer, and output layer. Input layer is a weight matrix designed to
pass the input data into the ANN through a linear mapping. In Hidden layers,
the input information is transformed through a matrix multiplication and an
activation function: $h_{1}=\sigma(W_{1}x+b_{1})$,
$h_{i}=\sigma(W_{i}h_{i-1}+b_{i}),i=\\{2,3,4,...,k\\}$, where
$x\in\mathbb{R}^{n_{x}}$ is the vector of input data,
$h_{i}\in\mathbb{R}^{n_{h_{i}}}$ is the $i_{th}$ hidden layer representation,
$W_{i}\in\mathbb{R}^{n_{h_{i}}\times{n_{h_{i-1}}}}$ and
$b_{i}\in\mathbb{R}^{n_{h_{i}}}$ are weights and bias for connecting the i-th
and (i-1)-th hidden representations, $k$ is the number of hidden layers, and
$\sigma$ is a nonlinear activation function that introduces nonlinearity into
the neural network. Specifically, ReLU activation function,
$\sigma(x)=max(0,x)$, is used in this paper. The last fully-connected hidden
layer does not impose the activation function. Then, a softmax layer
calculates the score/probability of each class via
$\frac{e^{f_{y_{i}}}}{\sum_{j}e^{f_{j}}}$, where $f_{y_{i}}$ is the score of
the correct category, and $f_{j}$ is the score of the $j_{th}$ catagory.
Finally, the loss is calculated based on the true label and regularization.
The objective of the network is to minimize the loss and maximum the accuracy
of both training and testing datasets.

Figure 1: ANN Structure

##  III Deep Compression

Neural networks usually require enormous memory to fund their massive
parameters. In this section, three deep compression techniques are discussed
to reduce the model size and keep high accuracy.

###  III-A Pruning

Network pruning utilizes the connectivity between neuron s. Small weights
under a predefined threshold are removed to reduce the model size. As shown in
Figure 2, the left figure is the original weights matrix. After pruning with a
threshold 0.1, the pruned weights are illustrated in the right figure. The
sparsity of the matrix reduces model size, and can speeds up the calculation
with hardware optimizations.

Figure 2: The left figure is the original weight matrix, and the right figure
is the matrix after pruning with a threshold 0.1. Highlighted weights are
removed.

###  III-B Clustering

Clustering is a technique that groups weights with close values. First,
weights are divided into a fixed number of groups. Second, values in the same
group will be reset as the same value and thus all weights in the same group
can be referenced to the same address to reduce model size. In Figure 3, the
left top figure is a weight matrix before clustering. Weights in this original
matrix are grouped into 6 groups, and all weights in the same group are each
assigned with a same value. The gradient of each weight is calculated through
back-propagation and the gradient of each group is the sum over all elements.
Finally, weights are updated through the cluster groups. Clustering
essentially minimizes the amount of references needed for a model, and thus
reduces the storage requirement of this model.

Figure 3: An illustration of the gradient update for clustered weights: the
top left figure is the original weight matrix colored according to clustered
groups; the top middle figure is the corresponding group indices; the bottom
left figure is the gradient of the weight matrix; the bottom middle figure is
the grouped gradients according to the centroids; the bottom right figure
shows the sum of gradients of each group; and the top right figure represent
the final gradient update of clusters.

###  III-C Quantization

In ANN, weights are usually stored a in high precision format such as float64
or float32, which consumes a relatively large space. In the quantization
process, the precision used to store weights is degraded in the exchange for a
lower space consumption. In Figure 4, weights in the left figure are stored in
one decimal, while the right figure reduces the storage requirement via
rounding off the decimal value.

Figure 4: The left figure is the weight matrix before quantization; the right
figure is the weight matrix after quantization.

##  IV Fault Detection for the Tennessee Eastman Chemical Process

The Tennessee Eastman process is a benchmark for fault detection [5]. Based on
a baseline ANN structure, various combinations of deep compression techniques
are applied to reduce the model size and achieve a high accuracy.

###  IV-A Tennessee Eastman Process Description

As summarized in Table I, the Tennessee Eastman process provides 52
measurements and 21 different fault types (including ”no fault” as faultNumber
0) in its dataset. Moreover, Chiang [6] and Zhang [7] pointed out that due to
the absence of observable change in the mean, variance and the higher order
variances in the data, we do not have enough information to detect Fault 3, 9,
and 15. Therefore, this paper follows the common rule to exclude the three
faults above for the consideration, resulting in 18 classification types
(including ”no failure” as faultNumber 0).

Fault Number | Process Variable | Type
---|---|---
faultNumber 1 | A/C feed ratio, B composition con- stant (stream 4) | Step
faultNumber 2 | B composition, A/C ratio constant (stream 4) | Step
faultNumber 3 | D feed temperature (stream 2) | Step
faultNumber 4 | Reactor cooling water inlet tem- perature | Step
faultNumber 5 | Condenser cooling water inlet tem- perature | Step
faultNumber 6 | A feed loss (stream 1) | Step
faultNumber 7 | C header pressure loss - reduced availability (stream 4) | Step
faultNumber 8 | A, B, C feed composition (stream 4) | Random variation
faultNumber 9 | D feed temperature (stream 2) | Random variation
faultNumber 10 | C feed temperature (stream 4) | Random variation
faultNumber 11 | Reactor cooling water inlet tem- perature | Random variation
faultNumber 12 | Condenser cooling water inlet tem- perature | Random variation
faultNumber 13 | Reaction kinetics | Slow drift
faultNumber 14 | Reactor cooling water valve | Sticking
faultNumber 15 | Condenser cooling water valve | Sticking
faultNumber 16 | Unknown | Unkown
faultNumber 17 | Unknown | Unkown
faultNumber 18 | Unknown | Unkown
faultNumber 19 | Unknown | Unkown
faultNumber 20 | Unknown | Unkown
TABLE I: Fault Types in the TE Process

###  IV-B Baseline Artificial Neural Network Models

In Tennessee Eastman process measurements, one entry in the dataset might
correspond to multiple fault types. Following the state-of-art method in [5],
we construct and train 18 different ANNs to detect each of the 18 fault types.
In baseline, each ANN has 6 layers in the structure of
$52-64-256-128-256-128-64-2$. The input layer dimension is 52, which
corresponds to the 52 measurement. The output layer has a dimension of 2,
which corresponds to the 2 possible outcomes (fault or no fault). The
resulting 18 models are served as the baseline for comparison with compressed
models.

###  IV-C Deep Compression Results

In this section, all 7 different combinations of pruning, clustering, and
quantization are applied to the 18 ANNs. First, each of the three compression
techniques is applied individually to exam the performance. Then, each two of
the compression techniques applied to the original 18 ANNs. At Last, all three
techniques are applied to the original ANNs to compare the compression rate
and accuracy change.

The compression results of ANNs obtained by applying each of the three
techniques individually are shown in Table II, III, and IV. As shown in Table
II, when only pruning is applied, the average compressed rate is 64.0% , and
the average accuracy change is -0.5% for all 18 different types of fault
labels. This result is promising considering the average accuracy only changed
slightly, while the average model size is significantly reduced by 64.0%.
Table III shows the compression result after clustering only, which achieves
76.0% average compressed rate and only -1.2% average accuracy change. This
method is the best among combinations that only apply a single compression
technique. Table IV presents the result after only quantization, which
achieves an average compression rate of 72.0%, and an average accuracy change
of -7.8%.

Fault |  ANN Size (bytes) | ANN Acc (%) |  Compressed Size (bytes) |  Compressed Acc (%) |  Compressed Rate (%) | Acc Change (%)
---|---|---|---|---|---|---
0 | 486924 | 94.4 | 175207 | 94.4 | 64.0 | 0.0
1 | 487120 | 99.4 | 175235 | 95.3 | 64.0 | -4.1
2 | 486853 | 99.6 | 174988 | 99.6 | 64.1 | 0.0
3 | 486931 | 96.1 | 175159 | 94.4 | 64.0 | -1.7
4 | 486901 | 95.2 | 175354 | 95.3 | 64.0 | 0.1
5 | 486753 | 95.6 | 174983 | 99.8 | 64.1 | 4.2
6 | 486914 | 99.8 | 174765 | 99.8 | 64.1 | 0.0
7 | 486992 | 97.9 | 175556 | 98.5 | 64.0 | 6.0
8 | 486789 | 95.3 | 175412 | 95.3 | 64.0 | 0.0
9 | 486782 | 94.4 | 175086 | 94.4 | 64.0 | 0.0
10 | 486737 | 95.3 | 175564 | 95.3 | 64.0 | 0.0
11 | 487124 | 97.6 | 175382 | 97.3 | 64.0 | -0.3
12 | 486838 | 96.4 | 175283 | 99.1 | 64.0 | 2.7
13 | 486819 | 97.1 | 175111 | 95.3 | 64.0 | -1.8
14 | 486844 | 96.8 | 175266 | 96.8 | 64.0 | 0.0
15 | 486899 | 96.4 | 175400 | 92.0 | 64.0 | -4.4
16 | 486816 | 95.3 | 175605 | 94.4 | 64.0 | -0.9
17 | 486844 | 97.6 | 175296 | 94.4 | 64.0 | -3.2
Average |  |  |  |  | 64.0 | -0.5

TABLE II: Compression results with pruning

Fault |  ANN Size (bytes) | ANN Acc (%) |  Compressed Size (bytes) |  Compressed Acc (%) |  Compressed Rate (%) | Acc Change (%)
---|---|---|---|---|---|---
0 | 486924 | 94.4 | 117313 | 94.4 | 75.9 | 0.0
1 | 487120 | 99.4 | 117725 | 99.1 | 75.8 | -0.3
2 | 486853 | 99.6 | 116627 | 99.4 | 76.0 | -0.2
3 | 486931 | 96.1 | 118057 | 94.4 | 75.8 | -1.7
4 | 486901 | 95.2 | 116922 | 94.4 | 76.0 | -0.8
5 | 486753 | 95.6 | 115861 | 99.6 | 76.2 | 4.0
6 | 486914 | 99.8 | 114344 | 94.4 | 76.5 | -5.4
7 | 486992 | 97.9 | 115685 | 94.9 | 76.2 | -3.0
8 | 486789 | 95.3 | 116442 | 94.4 | 76.1 | -0.9
9 | 486782 | 94.4 | 117599 | 94.4 | 75.8 | 0.0
10 | 486737 | 95.3 | 117798 | 94.4 | 75.8 | -0.9
11 | 487124 | 97.6 | 115934 | 94.7 | 76.2 | -2.9
12 | 486838 | 96.4 | 117672 | 94.4 | 75.8 | -2.0
13 | 486819 | 97.1 | 116671 | 94.4 | 76.0 | -2.7
14 | 486844 | 96.8 | 118284 | 96.3 | 75.7 | -0.5
15 | 486899 | 96.4 | 116166 | 94.7 | 76.1 | -1.7
16 | 486816 | 95.3 | 117600 | 94.4 | 75.8 | -0.9
17 | 486844 | 97.6 | 117014 | 94.5 | 76.0 | -3.1
Average |  |  |  |  | 76.0 | -1.2

TABLE III: Compression results with clustering

Fault |  ANN Size (bytes) | ANN Acc (%) |  Compressed Size (bytes) |  Compressed Acc (%) |  Compressed Rate (%) | Acc Change (%)
---|---|---|---|---|---|---
0 | 486924 | 94.4 | 137620 | 94.4 | 71.7 | 0.0
1 | 487120 | 99.4 | 137030 | 95.9 | 71.9 | -3.5
2 | 486853 | 99.6 | 137482 | 98.8 | 71.8 | -0.8
3 | 486931 | 96.1 | 136678 | 94.4 | 71.9 | -1.7
4 | 486901 | 95.2 | 137766 | 94.4 | 71.7 | -0.8
5 | 486753 | 95.6 | 137140 | 99.5 | 71.8 | 3.9
6 | 486914 | 99.8 | 136784 | 13.7 | 71.9 | -86.1
7 | 486992 | 97.9 | 137008 | 94.8 | 71.9 | -3.1
8 | 486789 | 95.3 | 137479 | 94.4 | 71.8 | -0.9
9 | 486782 | 94.4 | 137911 | 94.4 | 71.7 | 0.0
10 | 486737 | 95.3 | 137582 | 94.4 | 71.7 | -0.9
11 | 487124 | 97.6 | 137075 | 95.5 | 71.9 | -2.1
12 | 486838 | 96.4 | 137428 | 71.8 | 75.8 | -25.6
13 | 486819 | 97.1 | 137064 | 82.9 | 71.8 | -14.2
14 | 486844 | 96.8 | 137421 | 95.8 | 71.8 | -1.0
15 | 486899 | 96.4 | 137109 | 96.1 | 71.8 | -0.3
16 | 486816 | 95.3 | 137560 | 94.4 | 71.7 | -0.9
17 | 486844 | 97.6 | 137238 | 95.0 | 71.8 | -2.6
Average |  |  |  |  | 72.0 | -7.8

TABLE IV: Compression results with Quantization

Next, a comprehensive study of applying two compression techniques are
conducted. Table V shows that when utilizing both clustering and pruning, the
average compressed rate increases to 87.3%, and the average accuracy change is
only -1.9%. Moreover, the compression result is consistent among all fault
types. The variance of compressed rate and accuracy change are respectively
0.16 and 1.72, which shows that there are few fluctuations in compression rate
and accuracy change, thus the average compression rate and average accuracy
change are reliably representative. Table VI shows the compression result with
both pruning and quantization. The average compression rate is 88.1%, which is
about the same as the result in Table V. But the average accuracy change drops
to -5.3%, which is the largest drop among methods that apply two compression
techniques. Table VII presents the compression result with both clustering and
pruning. The average accuracy change is -1.8%, which is close to the average
accuracy change in Table V. The compression rate is 82.2%, which is 5.1%
smaller than the result in Table V. In general, applying two compression
techniques achieves better results than only applying one single technique.

Fault |  ANN Size (bytes) | ANN Acc (%) |  Compressed Size (bytes) |  Compressed Acc (%) |  Compressed Rate (%) | Acc Change (%)
---|---|---|---|---|---|---
0 | 486924 | 94.4 | 62014 | 94.4 | 87.3 | 0.0
1 | 487120 | 99.4 | 61788 | 95.3 | 87.3 | -4.1
2 | 486853 | 99.6 | 63234 | 95.3 | 87.0 | -4.3
3 | 486931 | 96.1 | 61694 | 94.4 | 87.3 | -1.7
4 | 486901 | 95.2 | 62078 | 95.3 | 87.3 | 0.1
5 | 486753 | 95.6 | 62995 | 95.3 | 87.0 | -0.3
6 | 486914 | 99.8 | 62065 | 95.3 | 87.3 | -4.5
7 | 486992 | 97.9 | 60908 | 95.3 | 87.5 | -2.6
8 | 486789 | 95.3 | 60772 | 95.3 | 87.5 | 0.0
9 | 486782 | 94.4 | 61664 | 94.4 | 87.3 | 0.0
10 | 486737 | 95.3 | 62178 | 95.3 | 87.2 | 0.0
11 | 487124 | 97.6 | 61671 | 95.3 | 87.3 | -2.3
12 | 486838 | 96.4 | 62292 | 95.3 | 87.2 | 0.0
13 | 486819 | 97.1 | 61002 | 95.3 | 87.5 | -1.8
14 | 486844 | 96.8 | 62443 | 92.1 | 87.2 | -4.7
15 | 486899 | 96.4 | 60988 | 94.4 | 87.5 | -2.0
16 | 486816 | 95.3 | 60709 | 94.4 | 87.5 | -0.9
17 | 486844 | 97.6 | 61086 | 94.4 | 87.5 | -3.2
Average |  |  |  |  | 87.3 | -1.9

TABLE V: Compression results with Pruning and Clustering

Fault |  ANN Size (bytes) | ANN Acc (%) |  Compressed Size (bytes) |  Compressed Acc (%) |  Compressed Rate (%) | Acc Change (%)
---|---|---|---|---|---|---
0 | 486924 | 94.4 | 57823 | 94.4 | 88.1 | 0.0
1 | 487120 | 99.4 | 58337 | 94.4 | 88.0 | -5.0
2 | 486853 | 99.6 | 57901 | 94.2 | 88.1 | -5.4
3 | 486931 | 96.1 | 58162 | 94.4 | 88.1 | -1.7
4 | 486901 | 95.2 | 58014 | 94.4 | 88.1 | -0.8
5 | 486753 | 95.6 | 57901 | 94.4 | 88.1 | -1.2
6 | 486914 | 99.8 | 58132 | 94.4 | 88.1 | -5.4
7 | 486992 | 97.9 | 57978 | 94.4 | 88.1 | -3.5
8 | 486789 | 95.3 | 57330 | 94.4 | 88.2 | -0.9
9 | 486782 | 94.4 | 58013 | 94.4 | 88.1 | 0.0
10 | 486737 | 95.3 | 58098 | 94.4 | 88.1 | -0.9
11 | 487124 | 97.6 | 57775 | 94.4 | 88.1 | -3.2
12 | 486838 | 96.4 | 58371 | 94.4 | 88.0 | -2.0
13 | 486819 | 97.1 | 56579 | 94.4 | 88.4 | -2.7
14 | 486844 | 96.8 | 58048 | 46.9 | 88.1 | -49.9
15 | 486899 | 96.4 | 57806 | 87.1 | 88.1 | -9.3
16 | 486816 | 95.3 | 57919 | 94.4 | 88.1 | -0.9
17 | 486844 | 97.6 | 57648 | 94.4 | 88.2 | -3.2
Average |  |  |  |  | 88.1 | -5.3

TABLE VI: Compression results with Pruning and Quantization

Fault |  ANN Size (bytes) | ANN Acc (%) |  Compressed Size (bytes) |  Compressed Acc (%) |  Compressed Rate (%) | Acc Change (%)
---|---|---|---|---|---|---
0 | 486924 | 94.4 | 86997 | 94.4 | 82.1 | 0.0
1 | 487120 | 99.4 | 86233 | 96.6 | 82.3 | -2.8
2 | 486853 | 99.6 | 86674 | 95.7 | 82.3 | -3.9
3 | 486931 | 96.1 | 85853 | 94.4 | 82.4 | -1.7
4 | 486901 | 95.2 | 86820 | 94.4 | 82.2 | -0.8
5 | 486753 | 95.6 | 86477 | 99.2 | 82.2 | 3.6
6 | 486914 | 99.8 | 85627 | 94.4 | 82.4 | -5.4
7 | 486992 | 97.9 | 86582 | 94.5 | 82.2 | -3.4
8 | 486789 | 95.3 | 86966 | 94.4 | 82.1 | -0.9
9 | 486782 | 94.4 | 86958 | 94.4 | 82.1 | 0.0
10 | 486737 | 95.3 | 86962 | 94.3 | 82.1 | -1.0
11 | 487124 | 97.6 | 86509 | 93.7 | 82.2 | -3.9
12 | 486838 | 96.4 | 86261 | 94.4 | 82.3 | -2.0
13 | 486819 | 97.1 | 86494 | 94.4 | 82.2 | -2.7
14 | 486844 | 96.8 | 86644 | 94.7 | 82.2 | -2.1
15 | 486899 | 96.4 | 86561 | 94.3 | 82.2 | -2.1
16 | 486816 | 95.3 | 86990 | 94.4 | 82.1 | -0.9
17 | 486844 | 97.6 | 86580 | 94.4 | 82.2 | -3.2
Average |  |  |  |  | 82.2 | -1.8

TABLE VII: Compression results with Clustering and Quantization

Finally, we exam the performance of applying all three compression techniques.
Table VIII shows that the average compressed rate increases to 91.5% while the
average accuracy change remains relatively small at -1.8%. The variances of
the compressed rates and accuracy changes are only 0.14 and 1.83,
correspondingly, which shows the consistency of the performance among all
fault types. Compare to results of applying two techniques in Table V, Table
VI, and Table VII, this method achieves the best compressed rate with no extra
loss in accuracy, which indicates applying all three techniques is better than
only applying two techniques for fault diagnosis in the TE process. Compared
to the results of applying only one compression technique in Table II, III,
and IV, although average accuracy change decreases slightly, applying all
three techniques still achieves accuracies higher than 94% for all fault
types. At the same time, the compressed rate 91.5% is significantly higher
than 76%, the best result of applying only one technique.

The results of all 7 different combination of compression techniques are
summarized in Figure 5. It is clear that applying all three techniques
achieves the highest compressed rate, as well as maintains a high average
accuracy above 94%. Thus, it is ideal to utilize all three compression
techniques for the fault diagnosis in the TE process.

Fault |  ANN Size (bytes) | ANN Acc (%) |  Compressed Size (bytes) |  Compressed Acc (%) |  Compressed Rate (%) | Acc Change (%)
---|---|---|---|---|---|---
0 | 486924 | 94.4 | 41391 | 94.4 | 91.5 | 0.0
1 | 487120 | 99.4 | 41264 | 94.4 | 91.5 | -5.0
2 | 486853 | 99.6 | 42833 | 94.4 | 91.2 | -5.2
3 | 486931 | 96.1 | 41320 | 94.4 | 91.5 | -1.7
4 | 486901 | 95.2 | 42055 | 95.3 | 91.4 | 0.1
5 | 486753 | 95.6 | 42466 | 95.5 | 91.3 | -0.1
6 | 486914 | 99.8 | 41165 | 95.3 | 91.5 | -4.5
7 | 486992 | 97.9 | 40674 | 95.3 | 91.6 | -2.6
8 | 486789 | 95.3 | 40523 | 95.3 | 91.7 | 0.0
9 | 486782 | 94.4 | 41471 | 94.4 | 91.5 | 0.0
10 | 486737 | 95.3 | 41642 | 95.3 | 91.4 | 0.0
11 | 487124 | 97.6 | 41741 | 95.3 | 91.4 | -2.3
12 | 486838 | 96.4 | 42006 | 95.3 | 91.4 | -1.1
13 | 486819 | 97.1 | 40762 | 95.3 | 91.6 | -1.8
14 | 486844 | 96.8 | 42055 | 95.3 | 91.4 | -1.5
15 | 486899 | 96.4 | 40559 | 94.4 | 91.7 | -2.0
16 | 486816 | 95.3 | 40741 | 94.4 | 91.6 | -0.9
17 | 486844 | 97.6 | 40543 | 94.4 | 91.7 | -3.2
Average |  |  |  |  | 91.5 | -1.8

TABLE VIII: Compression Statistics with Pruning, Clustering, and Quantization
Figure 5: The left plot shows compressed rates of all fault types over 7
different combinations of compression techniques, where P = pruning, C =
clustering, and Q = quantization. The right plot shows the average accuracy
over all fault types with 7 different combinations of compression techniques.
Results of different compression methods are colored consistently in two
plots.

##  V Conclusion

This paper studies deep compression techniques for fault diagnosis on the
Tennessee Eastman process. In respond to the demand of fast online detection
for industrial processes, three compression techniques (pruning, clustering,
and quantization) are applied to reduce model size and computational
complexity. We have examined a comprehensive list of 7 different combinations
of compression techniques. All methods achieve high model compression rates
over 64% while maintaining high fault detection accuracy. The best candidate
for fault detection on the Tennessee Eastman chemical process is applying all
three techniques, which reduces model sizes by 91.5% and remains a high
accuracy over 94%. This result leads to smaller storage requirement in
production environment, and makes the deployment smoother in real world.

## References

  * [1] Dan Ciresan, Ueli Meier, and Jürgen Schmidhuber.  Multi-column deep neural networks for image classification.  In IN PROCEEDINGS OF THE 25TH IEEE CONFERENCE ON COMPUTER VISION AND PATTERN RECOGNITION (CVPR 2012, pages 3642–3649, 2012.
  * [2] Richard Socher, Alex Perelygin, Jean Wu, Jason Chuang, Christopher D. Manning, Andrew Ng, and Christopher Potts.  Recursive deep models for semantic compositionality over a sentiment treebank.  In Proceedings of the 2013 Conference on Empirical Methods in Natural Language Processing, pages 1631–1642, Seattle, Washington, USA, October 2013. Association for Computational Linguistics.
  * [3] C. Chen, A. Seff, A. Kornhauser, and J. Xiao.  Deepdriving: Learning affordance for direct perception in autonomous driving.  In 2015 IEEE International Conference on Computer Vision (ICCV), pages 2722–2730, 2015.
  * [4] J.J. Downs and E.F. Vogel.  A plant-wide industrial process control problem.  Computers & Chemical Engineering, 17(3):245 – 255, 1993.  Industrial challenge problems in process control.
  * [5] Seongmin Heo and Jay H. Lee.  Fault detection and classification using artificial neural networks.  IFAC-PapersOnLine, 51(18):470 – 475, 2018.  10th IFAC Symposium on Advanced Control of Chemical Processes ADCHEM 2018\.
  * [6] Leo Chiang, E. Russell, and Richard Braatz.  Fault detection and diagnosis in industrial systems.  Measurement Science and Technology - MEAS SCI TECHNOL, 12, 10 2001\.
  * [7] Ying wei Zhang.  Enhanced statistical analysis of nonlinear processes using kpca, kica and svm.  Chemical Engineering Science, 64:801–811, 2009.