L PREPRINT - Accepted by 25th IEEE Latin American Test Symposium (LATS) 2024,
Maceio, Brazil3.5cm,1cm1.001.8

ReRAM

    Resistive Random Access Memory
CIM

    Computing-in-Memory
MVM

    Matrix-Vector Multiplication
NBB

    NeuroBreakoutBoard
NVM

    Non-Volatile Memories
TIA

    TransImpedance Amplifier
HRS

    High Resistive State
LRS

    Low Resistive State
WL

    Word Line
BL

    Bit Line
NN

    Neural Network
ADC

    Analog-to-Digital Converter
DAC

    Digital-to-Analog Converter
API

    Application Programming Interface
ISPVA

    Incremental Step Pulse Program and Verify Algorithm
NMOS

    Negative channel Metal-Oxide Semiconductor
MOS

    Metal-Oxide Semiconductor
CMOS

    Complementary Metal-Oxide Semiconductor
MIM

    Metal-Insulator-Metal
OxRAM

    Oxide-based Random Access Memory
SL

    Source Line
C2C

    Cycle-to-Cycle
D2D

    Device-to-Device
CDF

    Cumulative Distribution Function
CF

    Conductive Filament
WM

    Window Margin
IF

    Incremental Form
IFV

    Incremental Form and Verify
VCM

    Valence Change Memory
ECM

    Electrochemical Mechanism
TCM

    Thermochemical Mechanism
UART

    Universal Asynchronous Receiver/Transmitter
Protobuf

    Protocol Buffer
FMC

    FPGA Mezzanine Card

# A Fully Automated Platform for Evaluating ReRAM Crossbars

Rebecca Pelke1, Felix Staudigl1, Niklas Thomas1, Nils Bosbach1, Mohammed
Hossein1, Jose Cubero-Cascante1,
Leticia Bolzani Poehls2, Rainer Leupers1, and Jan Moritz Joseph1 This work was
partially funded by the Federal Ministry of Education and Research (BMBF,
Germany) in the project NeuroSys (Project Nos. 03ZU1106CA) and NEUROTEC
(project Nos. 16ES1134 and 16ES1133K).  1Institute for Communication
Technologies and Embedded Systems, RWTH Aachen University
2Chair of Integrated Digital Systems and Circuit Design, RWTH Aachen
University
{pelke, staudigl, thomas, bosbach, hossein, cubero, leupers<EMAIL_ADDRESS>aachen.de
<EMAIL_ADDRESS>
###### Abstract

rram is a promising candidate for implementing Computing-in-Memory (CIM)
architectures and neuromorphic circuits. Resistive Random Access Memory
(ReRAM) cells exhibit significant variability across different memristive
devices and cycles, necessitating further improvements in the areas of
devices, algorithms, and applications. To achieve this, understanding the
stochastic behavior of the different ReRAM technologies is essential. The
NeuroBreakoutBoard (NBB) is a versatile instrumentation platform to
characterize Non-Volatile Memories (NVMs). However, the NBB itself does not
provide any functionality in the form of software or a controller. In this
paper, we present a control board for the NBB able to perform reliability
assessments of 1T1R ReRAM crossbars. In more detail, an interface that allows
a host PC to communicate with the NBB via the new control board is
implemented. In a case study, we analyze the Cycle-to-Cycle (C2C) variation
and read disturb TiN/Ti/HfO2/TiN cells for different read voltages to gain an
understanding of their operational behavior.

###### Index Terms:

ReRAM, Test Automation, ReRAM variability

††publicationid: pubid: 979-8-3503-6555-9/24/$31.00 ©2024 European Union

##  I Introduction

There is an ever-growing demand for efficient Neural Network (NN) inference.
Due to the computationally and data-intensive workloads, traditional computers
quickly encounter the von Neumann bottleneck [1] and the memory-wall [2]. One
possible solution for addressing this problem is the adoption of CIM
architectures. The concept involves statically storing the weights of NNs in
non-volatile CIM cells. The computations take place within the memory,
tackling the challenges of traditional computing architectures [3, 4]. ReRAM
is a promising candidate for CIM due to its high device density, scalability,
switching performance, low power consumption, and Complementary Metal-Oxide
Semiconductor (CMOS) manufactuirng process compatibility [5, 6, 7].

The downside of ReRAM is poor uniformity between cycles and devices, known as
C2C and Device-to-Device (D2D) variabilities [8]. C2C variation is associated
with the stability and reliability of ReRAM devices over repeated cycles of
use. It is an intrinsic property of the technology [9, 10]. D2D variation is
associated with manufacturing and production consistency. Manufacturing
deviations in thickness, length, width, and surface roughness of the ReRAM
cells increase the D2D variation [11]. Note that these problems prevent the
rapid commercialization and adoption of ReRAM for implementing emerging
applications. To advance ReRAM technology, there are several topics to
explore, including

  1. 1.

enhancing ReRAM devices [12, 13, 14].

  2. 2.

improving program and read algorithms [15, 16, 17].

  3. 3.

increasing the fault tolerance of applications [18, 19, 20].

Regardless of the topic, understanding the characteristics of individual
devices is essential for properly implementing and optimizing these emerging
applications. Thus, automated measurement processes are required to reduce the
need for manual and time-consuming laboratory measurements. Different systems
have been developed for this purpose [21, 22, 23, 24], such as the NBB. The
NBB can measure a wider range of devices due to its flexibility [24].

The flexibility is achieved through a combination of the interconnect matrix,
custom-designed TransImpedance Amplifiers, and the omission of already-fixed
control mechanisms on the NBB. This means that the components of the NBB need
to be controlled by an external device because the NBB itself is not yet a
fully functional measurement system. We aim to close this gap by developing a
new platform able to perform reliability assessments of 1T1R ReRAM crossbars.
In this context, the main contributions of this paper are:

  * •

Development of a controller board that allows controlling the NBB from a host
PC, eliminating the need for supervision during measurements (increase of
measurement consistency and time savings in the laboratory).

  * •

Implementation of a new interface accessible from the host PC for easily
configuring the target experiments.

  * •

Validation of the proposed controller board and interface based on a case
study aiming the analysis of read voltage impact on C2C variability and read
disturb of TiN/Ti/HfO2/TiN cells.

The paper is organized as follows: Section II provides the required background
related to ReRAMs’ behavior as well as summarizes the main metrics adopted for
evaluating ReRAM cells. Section III describes the proposed controller board
and interface. Section IV summarizes the results related to the performed case
study’s reliability assessment. Finally, Section V draws a conclusion.

##  II Structure

This Section presents the background related to memristive devices. It further
introduces their main evaluation metrics.

### Memristive Devices: Structure and Functionality

A memristive device is composed of a transition-metal-oxide layer, e.g., HfO2,
ZrO2, Ta2O5, or TiO2, between two conducting electrodes [5]. After
manufacturing, especially in oxide-based filamentary-type devices, memristive
devices usually have very high electrical resistance. A large voltage is
required for the very first set operation, also known as the forming process
[25]. This process, a controlled soft breakdown, drastically reduces the
device resistance enabling the resistance-switching behavior in the subsequent
cycles for the filamentary regime [26]. The forming process directly impacts
the lifespan, forming yield, and C2C variability of the memristor. After
forming, set and reset operations can be applied, which bring the device to
the Low Resistive State (LRS) and High Resistive State (HRS), respectively
[27].

The authors of [28] compared three forming algorithms: Single pulse,
Incremental Form (IF), and Incremental Form and Verify (IFV). The first method
applies a single voltage pulse for forming memristors, the second one applies
multiple voltage pulses with increasing magnitudes, and the third, IFV,
involves increasing voltage pulses combined with verification steps to monitor
the state of the cell and to prevent current overshoots. They have found out
that IFV outperforms the other techniques in terms of cell yield and switching
behavior due to the tight control of the conductive filament creation.

In this work, the Incremental Step Pulse Program and Verify Algorithm (ISPVA)
[29], an IFV algorithm, is used. The algorithm, illustrated in Fig. 1, applies
voltage pulses with increasing magnitude to the electrodes to change the
resistive state. After each programming pulse $\mathrm{V_{prog}}$, a
verification voltage $\mathrm{V_{verify}}$ is applied to measure if the
desired current has been reached. The algorithm contains adjustable
parameters, e.g., the pulse duration $\mathrm{T_{pulse}}$, the start voltage
$\mathrm{V_{start}}$, the voltage increase $\mathrm{V_{step}}$, and the stop
voltage $\mathrm{V_{stop}}$. During set, positive voltages are applied to the
Word Line (WL) and $0\text{\,}\mathrm{V}$ to the Bit Line (BL). During reset,
positive voltages are applied to the BL and $0\text{\,}\mathrm{V}$ to the WL.

Figure 1: ISPVA [29] voltage sequence

###  ReRAMs: Structure and Metrics

In 1T1R crossbars, each ReRAM cell consists of one transistor and one
memristive device. The transistor is used as a switch to disconnect the
memristor from the crossbar, which reduces sneak-path currents [30]. There are
three types of 1T1R crossbars, namely typical, vertical, and pseudo crossbars
[31]. They differ in the horizontal or vertical placement of the Source Line
(SL), WL, and BL. The conductance value of the memristor $\mathrm{G_{ij}}$
represents the weight. The read voltage applied to the WL represents the
input, and the output current at the BL corresponds to the result of the
Matrix-Vector Multiplication (MVM) between the weights and inputs. The SL is
connected to the gate of the transistor. Fig. 2 shows a pseudo crossbar, which
is used later in the case study. In pseudo crossbars, WL and SL are both
placed in the horizontal direction.

Figure 2: 1T1R crossbar (pseudo crossbar)

The following properties are typically analyzed to characterize a ReRAM cell:

  * •

wm: $\mathrm{R_{OFF}/R_{ON}}$ ratio.

  * •

Retention: Stability of the cell (at higher temperatures).

  * •

Endurance: Maximum number of programming cycles $\mathrm{N_{cmax}}$ at a
constant Window Margin (WM).

  * •

Read disturb: Resistance-drifts after multiple readouts.

  * •

Variation: C2C variation (different cycles) and D2D variation (different
devices).

Different physical compositions of the ReRAM cells result in varying
behaviors. MOx/CuTex shows a good retention, while HfO2/CuTex reports large
WMs at a high endurance. HfO2/Ti present a good trade-off between endurance
and retention [32]. Related works have also revealed that programming
algorithms impact the switching behavior [33], variability, and reliability
[34, 16]. Further, the read voltages have significant effects on read disturb
and C2C variation. The authors of [17] analyzed this for Pt/ZrO2/Ta/Pt cells.
They have found out that it is advantageous to read in the reset direction. In
the following study, we plan to examine this phenomenon for HfO2/Ti cells. A
particular focus is placed on the HRS since a gradual HRS drift towards the
LRS is observed over cycling for this cell type [35]. When $\mathrm{N_{cmax}}$
is exceeded, programming the HRS becomes unfeasible, i.e., the cell is broken.

### The NeuroBreakoutBoard (NBB)

The aim of NBB is to provide a customizable infrastructure for measuring Non-
Volatile Memories (NVM) [24]. The NBB includes a configurable interconnect
matrix and circuitry for voltage pulse generation and current measurement.
These components are explained in more detail in the following.

####  II-1 Configurable Signal Interconnection Module

The board employs arrays of multiplexers to distribute different signal lines
or ground potential to each pin of the NVM socket independently. Around the
socket, $68$ $8$:$1$ multiplexers are spread over the four sides of the socket
(see Fig. 3). Each multiplexer receives $5$ signals from the Digital-to-Analog
Converter (DAC), a ground line, and an external voltage line. The south
multiplexers are additionally connected to the current sensing module.

####  II-2 Voltage Pulse Generator Module

A $12$-$\mathrm{bit}$,
$16$-$\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{n}\mathrm{e}\mathrm{l}$
DAC is used to generate the different voltage levels needed for the crossbar
operations. The $16$ channels can be used independently. Programming ReRAM
cells requires very short and precise voltage pulses [29]. For this purpose,
additional analog switches with switching times below $60\text{\,}\mathrm{ns}$
are provided.

####  II-3 Current Sensing Module

Accurate current sensing is required for different tasks like reading the
calculation result of a MVM. Furthermore, programming algorithms, such as
ISPVA, constantly require accurate feedback about the current. The currents
can be measured at the BLs of the crossbar (see Fig. 2). For accurate sensing,
the NBB has a circuit called TIA sense, which uses adjustable feedback
resistors to translate currents into voltages. The NBB has $16$ TIA sense
circuits, one for each south pin. To measure a wide range of currents,
different feedback resistances from $43\text{\,}\mathrm{\SIUnitSymbolOhm}$ to
$100\text{\,}\mathrm{k\SIUnitSymbolOhm}$ can be used. TIA sense further
includes an amplifier to adjust the voltage for the Analog-to-Digital
Converter (ADC). The ADC has $8$ $18$-$\mathrm{bit}$ channels, which can be
used simultaneously. The input resistance is
$1\text{\,}\mathrm{M\SIUnitSymbolOhm}$.

The authors of [24] evaluated the measurement accuracy of the NBB. For
measuring crossbar cells, we expect the resistances to range from
$1\text{\,}\mathrm{k\SIUnitSymbolOhm}$ to
$100\text{\,}\mathrm{k\SIUnitSymbolOhm}$. In this range, the measurement
accuracy remains highly reliable, with a relative error of less than
$0.5\text{\,}\mathrm{\char 37\relax}$.

##  III The Proposed Controler Board and Interface

Figure 3: Integration of NBB into automated measuring system

The main components of the controller board developed in this paper are shown
in Fig. 3. It can be plugged into the FPGA Mezzanine Card (FMC) connector of
the NBB. In general terms, the controller board generates the control signals
for the NBB modules. The STM32 microcontroller manages the configurable signal
interconnection module, the voltage pulse generator module, and the current
sensing module. In Fig. 3, the control parts are highlighted in the same color
as the corresponding modules of the NBB. All interfaces (host PC - controller
board - NBB) are marked in orange. The controller board enables the
communication between the host PC and the NVM chip. To be more precise, the
developed board

  * •

controls the multiplexers for each operation to connect the lines to the BLs,
WLs, and SLs of the NVM (blue).

  * •

implement ReRAM programming algorithms such as ISPVA by controlling the DAC
and analog switches. The voltage pulses must have different amplitudes,
shapes, and polarities depending on the operation and the underlying
technology (yellow).

  * •

control the feedback resistors and the TIA for accurate voltage sensing
(green).

### Host PC to Control Board Communication

For communication between the host PC and the control board, a Universal
Asynchronous Receiver/Transmitter (UART) interface with Protocol Buffer
(Protobuf) for message serialization is used. Before any command can be sent
to the control board, the setup phase of the controller must be completed to
initialize the SPI channels, the clock, and the GPIO interface. Then the NBB
is set up, i.e., the IO expanders are woken up, the DAC is cleared, the ADC
converter is initialized, the multiplexer arrays are reset by switching them
to ground potential, and the feedback resistors of the TIAs are set to their
highest values.

After the setup phase, the main loop starts on the microcontroller (see Fig.
4). Receive Request checks the host PC’s request for validity and unpacks the
message. If the operation is valid, it will be executed on the NVM chip in
Process Command. Finally, the result of the operation is sent back to the host
PC in Send Response. For example, this is the measured current of a read
operation.

Figure 4: Setup and main loop of the control board

### Python API

To abstract the Protobuf format and UART messages for the user, we developed a
host-side Application Programming Interface (API). It is a set of Python
functions that can be called to access the NVM cells. A sample of the API
functions is shown in Fig. 5. The corresponding cell can be selected with sl
and bl. Note that all API functions have a return value which is received by
Send Response (for the sake of clarity, this is not shown here). These Python
functions were also used to conduct the measurements presented in the case
study in Section IV.

⬇

iTargetSet=80, vGateRead=1.5)

resetCell(…, sl, bl, vGateReset=2.7,

iTargetReset=5, vGateRead=1.5)

setCell(…, sl, bl, vGate=1.5,

iTarget=80, vGateRead=1.5)

readCell(…, sl, bl, vGateRead=1.5,

vRead=vRead)

Figure 5: Extract of API functions

##  IV Validation of the Proposed Controller Board and Interface

((a)) Memristor 1: CDF of LRS and HRS, negative read
voltages$-140$$-120$$-100$$-80$$-60$$-40$$-20$$0$$0$$0.5$$1$Current
($\mathrm{\SIUnitSymbolMicro
A}$)CDF$-0.1\text{\,}\mathrm{V}$$-0.15\text{\,}\mathrm{V}$$-0.2\text{\,}\mathrm{V}$$-0.25\text{\,}\mathrm{V}$$-0.3\text{\,}\mathrm{V}$

((b)) Memristor 1: CDF of LRS and HRS, positive read
voltages$-10$$0$$10$$20$$30$$40$$50$$60$$70$$80$$90$$100$$110$$120$$130$$140$$0$$0.5$$1$Current
($\mathrm{\SIUnitSymbolMicro
A}$)CDF$0.1\text{\,}\mathrm{V}$$0.15\text{\,}\mathrm{V}$$0.2\text{\,}\mathrm{V}$$0.25\text{\,}\mathrm{V}$$0.3\text{\,}\mathrm{V}$

((c)) Memristor 2: CDF of LRS and HRS, negative read
voltages$-140$$-120$$-100$$-80$$-60$$-40$$-20$$0$$0$$0.5$$1$Current
($\mathrm{\SIUnitSymbolMicro A}$)CDF

((d)) Memristor 2: CDF of LRS and HRS, positive read
voltages$0$$20$$40$$60$$80$$100$$120$$140$$0$$0.5$$1$Current
($\mathrm{\SIUnitSymbolMicro A}$)CDF

Figure 6: CDF of LRS and HRS states for different cells using negative read
voltages (left) vs. positive read voltages (right).

To properly validate the proposed controller board and interface, we conduct a
case study using a $12\times 7$ TiN/Ti/HfO2/TiN pseudo crossbar. We evaluate
the C2C variation, read disturb, and endurance as metrics (see Section II).
The crossbar is plugged into the NBB’s socket. The ReRAM devices are formed
using the ISPVA [29] and the parameters provided by the vendor. After applying
a forming-reset-set sequence to all cells, the C2C variation and the read
disturb of the LRS and HRS are analyzed. In the following, the ISPVA algorithm
is also used to program the LRS and HRS (see Section II). Table I lists all
relevant parameters.

TABLE I: ISPVA set and reset parameters Parameter | Value
---|---
Top-electrode voltage LRS ($\mathrm{V_{te}}$) |  $0.5\text{\,}\mathrm{V}$ : $0.1$ |  $\mathrm{V}$ : $2\text{\,}\mathrm{V}$
Gate voltage LRS ($\mathrm{V_{g,lrs}}$) |  $1.5$ |  $\mathrm{V}$
Current target LRS ($\mathrm{I_{trg,lrs}}$) |  $80.0$ |  $\mathrm{\SIUnitSymbolMicro A}$
Bottom-electrode voltage HRS ($\mathrm{V_{be}}$) |  $0.5\text{\,}\mathrm{V}$ : $0.1$ |  $\mathrm{V}$ : $2\text{\,}\mathrm{V}$
Gate voltage HRS ($\mathrm{V_{g,hrs}}$) |  $2.7$ |  $\mathrm{V}$
Current target HRS ($\mathrm{I_{trg,hrs}}$) |  $5.0$ |  $\mathrm{\SIUnitSymbolMicro A}$

###  IV-A Cycle-To-Cycle Variation

We analyze the C2C variation for different read voltages, i.e., WL voltages
(see Fig. 2). We apply negative and positive read voltages, specifically
$-0.3\text{\,}\mathrm{V}$, $-0.25\text{\,}\mathrm{V}$,
$-0.2\text{\,}\mathrm{V}$, $-0.15\text{\,}\mathrm{V}$,
$-0.1\text{\,}\mathrm{V}$, $0.1\text{\,}\mathrm{V}$,
$0.15\text{\,}\mathrm{V}$, $0.2\text{\,}\mathrm{V}$,
$0.25\text{\,}\mathrm{V}$, and $0.3\text{\,}\mathrm{V}$. For each WL voltage,
we apply the set-read-reset-read sequence
$100\text{\,}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{s}$ and measure
the SL current using the measurement circuitry of the NBB (see Section III).

The results for two different cells are shown in Fig. 6. The Cumulative
Distribution Function (CDF) of the LRS and HRS for the same WL voltage share
the same color. In all measurements, both LRS and HRS show a higher C2C
variance at higher WL voltages (in terms of absolute value). It can be seen
that the measured currents in the HRS have a smaller variance for negative WL
voltages (Figs. 6(a) and 6(c)) than for positive ones (Figs. 6(b) and 6(d)).
The higher the (positive) WL voltage, the more likely the HRS swings towards
LRS when reading. This is because the read operation with positive voltages
has the same polarity as the set operation (see Table I).

The same behavior has been described by the authors of [17]. They suggest
reading in the reset direction because it results in a slower transition of
the resistance state, as opposed to reading in the set direction.
Additionally, employing lower voltages enhances the resistance stability.

###  IV-B Read Disturb

((a)) Read stability for cells with stable HRS

((b)) Read stability for cells with unstable HRS

Figure 7: Read stability for different ReRAM cells using different read
voltages

Since a major application of ReRAM crossbars is the inference of NNs, it is
important to evaluate the read disturb, which refers to the alteration of data
caused by read operations. A major advantage of CIM is that the weights,
represented as LRS and HRS, are programmed only once and then used for
multiple inferences. Therefore, the cells must remain stable even when read
several times. We analyze the read disturb for the WL voltages mentioned in
Section IV-A. For each of the $10$ WL voltages, we apply the sequence set-
read(100)-reset-read(100),
$50\text{\,}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{s}$, i.e, $1,000$
write operations and $100,000$ read operations are performed per cell.

Fig. 7 shows the read stability of a stable and an unstable cell. Positive BL
currents occur as a result of positive WL voltages and negative currents as
caused by negative WL voltages. For clarity, only every $10$th read cycle is
shown. The box plots illustrate the variation of $50$ independent
measurements. The box boundaries represent the $25\text{\,}\mathrm{\char
37\relax}$ and $75\text{\,}\mathrm{\char 37\relax}$ quartiles. The lower and
upper whiskers denote the $2.5\text{\,}\mathrm{\char 37\relax}$ and
$97.5\text{\,}\mathrm{\char 37\relax}$ boundaries, respectively. The median is
highlighted in orange, and the dashed lines connect the mean values. The HRS
is always the state closer to the x-axis (see also Fig. 6). We show results
for cells with a stable (Fig. 7(a)) and unstable (Fig. 7(b)) HRS. For both
cases, negative read voltages, i.e., voltages in the reset direction, are
preferred. In cells with an unstable HRS, the cell is continuously
reprogrammed and switches to the LRS after just a single-digit read operation.
This effect becomes more pronounced with higher WL voltages. As seen in the
example from Fig. 7(b), HRS and LRS states significantly overlap with positive
WL voltages after only $10$ read accesses. For read operations with negative
WL voltages, all states remain distinguishable after $50$ accesses. Despite
separable states, the high variance in HRS can harm computing operations.

The maximum endurance of a particular ReRAM technology depends on the maximum
programming energy the system can endure before failure [32]. This means that
each set and reset operation contributes to the gradual degradation of the
oxide material, leading to failure at $N_{cmax}$. The reported maximum
endurance for Oxide-based Random Access Memory (OxRAM) technology has been
documented within the range of $10^{3}$ to $10^{9}$ cycles [32]. In the case
of HfO2/Ti cells, oxide degradation primarily impacts the HRS, resulting in a
gradual drift of the HRS towards the LRS [35]. This explains the behavior
observed in Fig. 7(b). The observed $N_{cmax}$ values of the tested crossbar
exhibited significant variability among different cells. While some cells have
a small $N_{cmax}$ of $100$ measurements, others show a stable HRS even after
$10,000$ cycles.

##  V Conclusion

In this work, we developed a fully automated test platform to characterize
ReRAM cells in 1T1R crossbar arrays. We implemented a control board that
enables communication between a host PC and the NBB. Using our Python-based
API, we are able to automate test cases for cell measurements.

The proposed controller board and interface were validated in a case study.
More specifically, we analyzed the effect of read voltages on C2C variation
and read disturb in TiN/Ti/HfO2/TiN cells. The findings reveal that reading in
reset direction achieves better results, especially for the read disturb, as
reading in set direction tends to swing the HRS towards the LRS.

Finally, this work optimizes laboratory measurement procedures, offering
considerable time savings. Cell measurements, which may span several hours,
can now be conducted without continuous supervision. Additionally, automated
measurements ensure the consistency of data collection by adhering to
predefined and standardized procedures, minimizing variations and human
errors. Looking ahead, the next steps are to add support for MVMs and to
automate test cases for multi-bit measurements.

## References

  * [1] X. Zou, S. Xu, X. Chen, L. Yan, and Y. Han, “Breaking the von Neumann bottleneck: architecture-level processing-in-memory technology,” _Science China Information Sciences_ , vol. 64, no. 6, p. 160404, 2021.
  * [2] X. Huang, C. Liu, Y.-G. Jiang, and P. Zhou, “In-memory computing to break the memory wall,” _Chinese Physics B_ , vol. 29, no. 7, p. 078504, 2020.
  * [3] S.-T. Wei, B. Gao, D. Wu, J.-S. Tang, H. Qian, and H.-Q. Wu, “Trends and challenges in the circuit and macro of RRAM-based computing-in-memory systems,” _Chip_ , vol. 1, no. 1, p. 100004, 2022.
  * [4] A. Baroni, A. Glukhov, E. Pérez, C. Wenger, D. Ielmini, and C. Zambelli, “An energy-efficient in-memory computing architecture for survival data analysis based on resistive switching memories,” _Frontiers in Neuroscience_ , 2022.
  * [5] J. Hazra, M. Liehr, K. Beckmann, M. Abedin, S. Rafq, and N. Cady, “Optimization of Switching Metrics for CMOS Integrated HfO2 based RRAM Devices on 300 mm Wafer Platform,” in _2021 IEEE International Memory Workshop (IMW)_. IEEE, 2021, pp. 1–4.
  * [6] J. S. Meena, S. M. Sze, U. Chand, and T.-Y. Tseng, “Overview of emerging nonvolatile memory technologies,” _Nanoscale research letters_ , 2014.
  * [7] J. S. Vetter and S. Mittal, “Opportunities for nonvolatile memory systems in extreme-scale high-performance computing,” _Computing in Science and Engineering_ , 2015.
  * [8] N. Raghavan _et al._ , “Stochastic variability of vacancy filament configuration in ultra-thin dielectric RRAM and its impact on OFF-state reliability,” in _2013 IEEE International Electron Devices Meeting_ , 2013\.
  * [9] H.-S. P. Wong _et al._ , “Metal–oxide RRAM,” _Proceedings of the IEEE_ , 2012.
  * [10] S. Yu, X. Guan, and H.-S. P. Wong, “On the stochastic nature of resistive switching in metal oxide RRAM: Physical modeling, Monte Carlo simulation, and experimental characterization,” in _2011 International Electron Devices Meeting_. IEEE, 2011.
  * [11] E. Brum _et al._ , “Evaluating the Impact of Process Variation on RRAMs,” in _LATS_. IEEE, 2021.
  * [12] Y. Fang _et al._ , “Improvement of HfO x-based RRAM device variation by inserting ALD TiN buffer layer,” _IEEE Electron Device Letters_ , 2018.
  * [13] W. Wu _et al._ , “A methodology to improve linearity of analog RRAM for neuromorphic computing,” in _2018 IEEE symposium on VLSI technology_. IEEE, 2018.
  * [14] T. Kempen, R. Waser, and V. Rana, “50x Endurance improvement in TaOx RRAM by extrinsic doping,” in _IMW_. IEEE, 2021.
  * [15] H. Liu _et al._ , “Uniformity improvement in 1T1R RRAM with gate voltage ramp programming,” _IEEE Electron Device Letters_ , 2014.
  * [16] A. Baroni, C. Zambelli, P. Olivo, E. Pérez, C. Wenger, and D. Ielmini, “Tackling the Low Conductance State Drift through Incremental Reset and Verify in RRAM arrays,” in _IIRW_. IEEE, 2021.
  * [17] C. Bengel _et al._ , “Reliability aspects of binary vector-matrix-multiplications using ReRAM devices,” _Neuromorphic computing and engineering_ , 2022.
  * [18] L. Xia, M. Liu, X. Ning, K. Chakrabarty, and Y. Wang, “Fault-tolerant training with on-line fault detection for RRAM-based neural computing systems,” in _DAC_ , 2017.
  * [19] L. Xia _et al._ , “Stuck-at fault tolerance in RRAM computing systems,” _IEEE JETCAS_ , 2017.
  * [20] L. Chen _et al._ , “Accelerator-friendly neural-network training: Learning variations and defects in RRAM crossbar,” in _DATE_. IEEE, 2017.
  * [21] R. Berdan, A. Serb, A. Khiat, A. Regoutz, C. Papavassiliou, and T. Prodromakis, “A $\mu$-Controller-Based System for Interfacing Selectorless RRAM Crossbar Arrays,” _IEEE Transactions on Electron Devices_ , 2015.
  * [22] J. Cayo and I. Vourkas, “Design Steps towards a MCU-based Instrumentation System for Memristor-based Crossbar Arrays,” in _MOCAST_. IEEE, 2021.
  * [23] R. De La Fuente, I. Vourkas, and M. Perez, “On the Development of MCU-based ad hoc HW Interface Circuitry for Memristor Characterization,” in _ECCTD_. IEEE, 2020.
  * [24] F. Staudigl _et al._ , “Work-in-Progress: A Universal Instrumentation Platform for Non-Volatile Memories,” in _ESWEEK_ , 2023.
  * [25] Guitarra, Silvana and Trojman, Lionel and Raymond, Laurent, “Resistive Switching Model of OxRAM Devices Based on Intrinsic Electrical Parameters,” in _2019 Latin American Electron Devices Conference (LAEDC)_ , 2019.
  * [26] L. B. Poehls _et al._ , “Review of Manufacturing Process Defects and Their Effects on Memristive Devices,” _Journal of electronic testing_ , 2021\.
  * [27] F. M. Puglisi, P. Pavan, A. Padovani, and L. Larcher, “A study on HfO2 RRAM in HRS based on I–V and RTN analysis,” _Solid-State Electronics_ , 2014\.
  * [28] A. Grossi _et al._ , “Electrical characterization and modeling of pulse-based forming techniques in RRAM arrays,” _Solid-State Electronics_ , vol. 115, pp. 17–25, 2016.
  * [29] E. Perez, C. Zambelli, M. K. Mahadevaiah, P. Olivo, and C. Wenger, “Toward Reliable Multi-Level Operation in RRAM Arrays: Improving Post-Algorithm Stability and Assessing Endurance/Data Retention,” _IEEE Journal of the Electron Devices Society_ , 2019.
  * [30] M. Mao, Y. Cao, S. Yu, and C. Chakrabarti, “Optimizing Latency, Energy, and Reliability of 1T1R ReRAM Through Cross-Layer Technique,” _IEEE Journal on Emerging and Selected Topics in Circuits and Systems_ , 2016.
  * [31] S. Son, C. La Torre, A. Kindsmüller, V. Rana, and S. Menzel, “A Study of the Electroforming Process in 1T1R Memory Arrays,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , 2022.
  * [32] C. Nail _et al._ , “Understanding RRAM endurance, retention and window margin trade-off using experimental results and simulations,” in _2016 IEEE International Electron Devices Meeting (IEDM)_. IEEE, 2016.
  * [33] E. Perez, M. K. Mahadevaiah, E. P.-B. Quesada, and C. Wenger, “Variability and energy consumption tradeoffs in multilevel programming of RRAM arrays,” _IEEE Transactions on Electron Devices_ , vol. 68, no. 6, pp. 2693–2698, 2021\.
  * [34] A. Baroni _et al._ , “Low Conductance State Drift Characterization and Mitigation in Resistive Switching Memories (RRAM) for Artificial Neural Networks,” _IEEE Transactions on Device and Materials Reliability_ , 2022\.
  * [35] G. Sassine _et al._ , “Optimizing Programming Energy for Improved RRAM Reliability for High Endurance Applications,” in _2018 IEEE International Memory Workshop (IMW)_. IEEE, 2018.

  *[CIM]: Computing-in-Memory
  *[ReRAM]: Resistive Random Access Memory
  *[NBB]: NeuroBreakoutBoard
  *[C2C]: Cycle-to-Cycle
  *[NN]: Neural Network
  *[NNs]: Neural Network
  *[CMOS]: Complementary Metal-Oxide Semiconductor
  *[D2D]: Device-to-Device
  *[LRS]: Low Resistive State
  *[HRS]: High Resistive State
  *[IF]: Incremental Form
  *[IFV]: Incremental Form and Verify
  *[ISPVA]: Incremental Step Pulse Program and Verify Algorithm
  *[WL]: Word Line
  *[BL]: Bit Line
  *[ReRAMs]: Resistive Random Access Memory
  *[SL]: Source Line
  *[MVM]: Matrix-Vector Multiplication
  *[WM]: Window Margin
  *[WMs]: Window Margin
  *[NVM]: Non-Volatile Memories
  *[DAC]: Digital-to-Analog Converter
  *[BLs]: Bit Line
  *[TIA]: TransImpedance Amplifier
  *[ADC]: Analog-to-Digital Converter
  *[FMC]: FPGA Mezzanine Card
  *[WLs]: Word Line
  *[SLs]: Source Line
  *[UART]: Universal Asynchronous Receiver/Transmitter
  *[Protobuf]: Protocol Buffer
  *[TIAs]: TransImpedance Amplifier
  *[API]: Application Programming Interface
  *[CDF]: Cumulative Distribution Function
  *[OxRAM]: Oxide-based Random Access Memory
  *[MVMs]: Matrix-Vector Multiplication


# Predicting malaria dynamics in Burundi using deep Learning Models

Daxelle Sakubu<EMAIL_ADDRESS>Kelly Joelle Gatore Sinigirira  David
Niyukuri  Doctoral School, University of Burundi, Burundi  Department of
Mathematics, University of Burundi, Burundi  The South African Department of
Science and Technology–National Research Foundation (DST-NRF) Centre of
Excellence in Epidemiological Modelling and Analysis (SACEMA), Stellenbosch
University, Cape Town, South Africa

###### Abstract

Malaria continues to be a major public health problem on the African
continent, particularly in Sub-Saharan Africa. Nonetheless, efforts are
ongoing, and significant progress has been made. In Burundi, malaria is among
the main public health concerns. In the literature, there are limited
prediction models for Burundi. We know that such tools are much needed for
interventions design. In our study, we built machine-learning based models to
estimates malaria cases in Burundi. The forecast of malaria cases was carried
out at province level and national scale as well. Long short term memory
(LSTM) model, a type of deep learning model has been used to achieve best
results using climate-change related factors such as temperature, rainfal, and
relative humidity, together with malaria historical data and human population.
With this model, the results showed that at country level different tuning of
parameters can be used in order to determine the minimum and maximum expected
malaria cases. The univariate version of that model (LSTM) which learns from
previous dynamics of malaria cases give more precise estimates at province-
level, but both models have same trends overall at provnce-level and country-
level.

###### keywords:

Malaria , Prediction , Deep Learning , Long Short Term Memory (LSTM) , Burundi

††journal: Epidemics

##  1 Introduction

Malaria is an infectious disease caused by the falciparum parasite, which is
transmitted through the bite of a female mosquito named Anopheles. Malaria
parasite are principally classified into four type: Plasmodium Falciparum,
Plasmodium Vivax, Plasmodium Malaria, and Plasmodium Ovale. Plasmodium
Falciparum is the most dangerous parasite since its incubation period is short
6-14 days[1] and it’s the principal cause of Malaria in Burundi nearly 90 %
[2].

According to the World Health Organization (WHO) report [3], estimated malaria
cases increased from 213 million to 228 million, and deaths increased from 534
000 to 602 000 in the WHO African Region between 2019 and 2020. This region
accounted for approximately 95 % of all cases and 96 % of all deaths globally,
children under the age of five account for 80 % of all deaths in this region.

In Burundi, malaria infection is among main major public health concerns after
Diarrheal diseases, Neonatal disorders, and tuberculosis [4, 5, 3, 6]. During
the recent malaria outbreak in 2017 [7, 8, 9], the World Vision International
assumed that climate change, population density, changed agricultural methods,
food hardship, and a lack of information and action to prevent malaria were
the main driving factors that contribute to the disease’s prevalence in
Burundi [10].

On a global scale, coordinated initiatives to minimize the malaria epidemic
are being planned as part of the millennium development goals. In Burundi, the
package of interventions for malaria control during pregnancy comprises ITN
promotion and use, IPTp with sulfadoxine-pyrimethamine under directly observed
treatment, and quick and successful treatment. The NMCP has yet to follow WHO
recommendations from 2016, which increased the recommended number of prenatal
care visits from four to eight. With Malaria Operational Plan (MOP) funds from
Fiscal Year (FY) 2021 [11], the team proposes to trial an evidence-based group
prenatal care strategy to enhance IPTp uptake (according to DHIS2, 54 % of
women received the recommended three doses at U.S. Government-supported health
centers in the first quarter of FY 2021). Although massive effort have already
been deployed, Malaria still prevail. Thus this project aim to understand the
climatic contributions by using machine learning models to predict Malaria
cases.

In recent years, there has been a great deal of interest in the development
and application of machine learning (ML) in the field of infectious diseases
[12]. Not only as a catalyst for academic studies, but also as a critical
means of detecting pathogenic microorganisms, implementing public health
surveillance, investigating host-pathogen interactions, discovering drug and
vaccine candidates, and so on. According to one survey, ML is used in 77% of
the products we use today. Machine learning (ML) is a subfield of Artificial
Intelligence (AI) that is an important tool in bioinformatics [13]. When
confronted with a range of large and complex data sets that must be analyzed,
ML may employ sophisticated algorithms and efficient models to extract
meaningful information from vast amounts of complex data-sets [14, 15, 16] .

Machine learning extracts useful information from enormous amounts of data by
using algorithms to recognize patterns and learn in an iterative process.
Instead of depending on any preconceived equation that may serve as a model,
ML algorithms use computing methods to learn directly from data [17]. The
union of mathematics and computer science in ML has shown significant
potential as a breakthrough in science and technology, and it has been applied
to a wide range of scientific fields, including biology. In the case of
Malaria prediction in Burundi some studies have already been done, and one has
used machine learning such as artificial neural network for prediction
[18].The authors investigated malaria across different group of ages and the
impact of meteorological factors on the high number of malaria cases during
some seasons.

Different degrees of accuracy were reported from previous investigations. The
overarching goal of this research is to investigate malaria cases predictions
based on meteorological data. Deep learning model will be used to forecast
malaria incidence in all five provinces. This study may help in the
improvement of public health measures especially on district level.

##  2 Data Description

The study was carried out using monthly data, collected from different
sources, namely: Geographical Institute of Burundi (Institut Géographique du
Burundi, IGEBU), the Institute of Statistics and Economic Studies of Burundi
Institut de Statistiques et d’Études Économiques du Burundi, ISTEEBU) and
Burundi National Malaria Control Program (NMCP). The data was collected for
all the eighteen provinces of the country.

###  2.1 Data extraction

Data collected from IGEBU was on a monthly scale from 2010 to 2022 with
parameters such as relative humidity, rainfall and temperature with their
maximum and minimum values. The average was calculated and inserted in the
data-set. Malaria historical data was obtained at NMCP Burundi on a monthly
scale from January 2010 to December 2022 for all previously eighteen
provinces. The human population feature was available online on ISTEEBU
website and the human population was calculated annually.

###  2.2 Data processing

The data collected was on different time scale monthly and annually, thus the
human population was considered as constant during the whole year. The
meteorological data contained some missing values that were filled using an
algorithm called missForest [19]. This algorithm was judged to be suitable
since it takes into account the possible relation between variables. Recently,
the communes, zones and hills/neighbourhoods of the Republic of Burundi have
been the subject of administrative redistricting-regrouping. Bujumbura
(Western part), Buhumuza (Eastern part), Gitega (Central part), Burunga
(southern part) and Butanyerera (Northern part) are the five provinces of the
country in the new delimitation of provinces.

In this reform, the provinces were regrouped as follow: Bujumbura included all
the commune of Bujumbura Mairie,Bujumbura Rural, Bubanza and Cibitoke. Gitega
assembled all the commune of Gitega, Mwaro, Karuzi and Muramvya. Buhumuza
grouped together all the communes of Cankuzo, Muyinga and Ruyigi. Butanyerera
aggregated all the communes of Kirundo, Ngozi and Kayanza. Finally, Burunga
included all the communes of Bururi, Makamba, Rumonge and Rutana. Therefore
,the data-set of these new delimited provinces were the mean of the old
provinces regrouped per month for the meteorological data and the sum for the
human population and malaria cases data.

##  3 Methodology

After processing the data, we start building the neural network models.The
prediction model which was build with four layers and several units.
Nevertheless, the performance were not good, hence a different approach was
taken in order to find good results.

###  3.1 Building ML models

Since the data-set were on a monthly scale, a different model that took into
account the previous information seems to be appropriate for this study. The
kind of algorithm that take into account previous information into a
chronological manner are recurrent neural network. This kind of neural network
in contrary to others have feedback connections. The network’s connection
weights and biases change once per training episode, similar to how
physiological changes in synaptic strengths store long-term memories; the
network’s activation patterns change once per time-step, similar to how the
brain’s electrical firing patterns change moment to moment to store short-term
memories. Recurrent neural networks are frequently used in multiple domain to
predict future events based on prior experience.

For instance, LSTM can be used for tasks like connected, unsegmented
handwriting recognition, video games,speech recognition,automated translation,
healthcare, Speech activity detection and Robotics. Since LSTM has proven to
have a good performance to aid with rational decision-making, it was utilized
for prediction of malaria cases in this study.

Figure 1:  LSTM cell

In Figure 1, a long short term memory cell is depicted where the square
represent the layers , the ellipse are the component wise operation,Ci are
hidden state vectors and Xi are input vector to the LSTM unit and hi are
hidden state vector also known as output vector (yi) of the LSTM unit.

Figure 2: Univariate LSTM.

In Figure 2, a univariate LSTM model is shown in a sequential manner with
previous states along side the input vectors $x_{i}$ and outputs $y_{i}$. The
previous malaria cases are the input vector and the output vector are the
actual malaria cases.

The Multivariate LSTM model is represented in Figure 3 with the multiple input
vector $x_{i}$ up to $x_{n}$ and the target $y_{i}$. In the multivariate LSTM
model the input vectors are the climate data human population and previous
malaria cases, the output vector are the current malaria cases.

Figure 3: Multivariate LSTM.

After tuning the data-set in the ML models we obtained results which are
presented in the following section.

##  4 Results

After running and fitting the data to the models, the error between the actual
and predicted cases were calculated using root mean square error(RMSE). From
the obtained results are observed in Table 1, overall,the univariate LSTM
model got the smallest RMSE while the Multivariate got the biggest RMSE. In
Bujumbura province the multivariate error is near four times as big as the
univariate error. In the case of Gitega province the error gap is three times
big between multivariate and univariate and this is the same for Butanyerera
province as well. In Burunga province the multivariate error is twice bigger
the univariate error. For Buhumuza province the error gap was huge, the
multivariate error was about eight times bigger than the univariate one. On
the country level prediction the multivariate error disparity was three times
big the univariate one.

Table 1: RMSE of malaria cases between Univariate LSTM and Multivariate LSTM. Province | Univariate LSTM | Multivariate LSTM
---|---|---
Bujumbura | $4868.69$ | $16777.17$
Gitega | $10943.18$ | $31012.00$
Burunga | $6403.33$ | $12964.48$
Butanyerera | $8288.07$ | $25187.09$
Buhumuza | $5664.18$ | $44893.25$
Country level: Burundi | $31635.95$ | $119724.68$

###  4.1 Province-level predictions

The experiment were done on eighty percent of the data-set while the testing
was done on the rest of the data-set. Malaria cases prediction on province
level are shown below in figure 5 from 2020 september to september 2022. In
univariate LSTM prediction, the curve trend are followed in most cases with
the observed cases being sligtly higher than the expected ones. In the
multivariate LSTM prediction, the curve trend are not coordinated except for
the country level prediction. For Gitega province the overall malaria cases
during that period of time are 2395666, the univariate model predicted
2444657,664 while the multivariate predicted 2338582,423. In Burunga province
the observed malaria cases were 1901832 while the univariate predicted were
1890808,86 and the multivariate 1974578,2. In Buhumuza province, the observed
cases during that time was 2772583 whereas the Univariate model predicted
2808850,717 and the multivariate 3567723,29. For Butanyerera, the actual cases
were 3568702 when the univariate predicted 3507692,91 and the multivariate
3808420,97. The Bujumbura province observed malaria cases were 2253588 during
that amount of time when the expected was 2232670,579 for Univariate LSTM
model and 2280516,361 for the multivariate LSTM model.

Figure 4: Predictions of malaria cases in provinces using Univariate LSTM
(left) and Multivariate LSTM (right)

###  4.2 Country-level predictions

After predicting malaria cases on province level, the country level prediction
was done. The results are depicted in Figure 5.The data-set of the country
level was the mean of all provinces for the climate data and the sum for the
human population and malaria cases. The curve trends were coordinated to the
actual trend in all models. Neverthless, the multivariate predicted more
malaria cases than observed one while the univariate predicted less. The
malaria cases that occurred in the country during the two years were
12959182,46, the univariate model predicted 12841653,9 which is around one
million cases less than what was observed. The multivariate LSTM model
predicted 15215766,15 that is around more than two millions and half than the
actual cases.

Figure 5: Predictions of malaria cases on country level using Univariate LSTM
(left) and Multivariate LSTM (right)

##  5 Discussion

Based on the climate related variables, the LSTM model provided promising
results on the prediction of malaria cases at country level. At province
level, for the north-eastern part Butanyerera and Buhumuza, the multivariate
LSTM tend to predict much more than observed. However, for Bujumbura and
Burunga, the south-Eastern part , the trend were followed up until october
2021 and start to decrease in 2022. In Gitega province the central part, the
expected cases from multivariate LSTM model were highly inferior to the
observed cases. The univariate LSTM on the other hand gave the best accuracy
with the curve trend being followed on country level as well as on the
province level.

This model proved to be useful in prediction of possible upcoming outbreaks.
However, few points need to be underligned. Since the the models use previous
information in order to predict the upcoming event, information gathering and
sharing is critical, without the most update information the accuracy will not
be the same and may decrease.

From Figure 5, we have seen that the Univariate LSTM results depict the
overall trends of malaria cases, and it gives more precise estimates at
province-level. With Univariate LSTM getting the best accuracy, multivariate
LSTM may give the maximum cases that can be expected overall. The models got
the same trends overall, especially on a country level prediction, this get to
confirm how much the climate condition affect the malaria outbreak in the
country.

The effects of climate change on disease outbreaks go beyond malaria
transmission dynamics. Investigating the meteorological factors in the case of
malaria alongside historical data such as human population and total malaria
cases per month was done in this work. With the results obtained, it’s clear
that despite artificial neural network acting like black box,there is a room
for meteorological data being an important factor to be considered in the
prediction of an outbreak especially on the country level where the maximum
cases were predicted. The use of deep learning model such as recurrent neural
networks has provided promising results in the prediction of malaria which
still causes outbreak in sub-saharan Africa. More applications of artificial
intelligence will help the collaborations between different intervention teams
in the healthcare system, especially on how they should respond after
predicting potential increase or decrease of malaria cases.

## References

  * [1] T Sajana and MR Narasingarao.  Machine learning techniques for malaria disease diagnosis-a review.  J Adv Res Dyn Control Syst, 9(6):349–369, 2017.
  * [2] Medecins Sans Frontieres.  More about malaria in Burundi, 2017.  Accessed: 2021-05-26.
  * [3] World Health Organization (WHO), Geneva.  World malaria report 2021, 2022.
  * [4] World Health Organization (WHO), Geneva.  World malaria report 2016, 2017.
  * [5] World Health Organization (WHO), Geneva.  World malaria report 2017, 2018.
  * [6] Health metrics for Burundi.  http://www.healthdata.org/burundi, 2020.  Accessed: 2021-05-26.
  * [7] World Health Organization (WHO)/Regional Office for Africa.  Weekly Bulletin on Outbreaks and other Emergencies: Week 27: 01–07 July 2017, 2017.
  * [8] United Nations For Children (UNICEF)/Burundi.  UNICEF Burundi Humanitarian Situation Report – 31 March 2017, 2017\.
  * [9] Pat Lok and Stijntje Dijk.  Malaria outbreak in burundi reaches epidemic levels with 5.7 million infected this year.  British Medical Journal Publishing Group, 2019.
  * [10] World Vision.  Eight facts about Burundi’s malaria epidemic, 2017.  Accessed: 2021-05-26.
  * [11] USAID.  BURUNDI MALARIA OPERATIONAL PLAN Fiscal Year 2021, 2022.  Accessed: 2022-05.
  * [12] Timothy L Wiemken and Robert R Kelley.  Machine learning in epidemiology and health outcomes research.  Annual review of public health, 41:21–36, 2019.
  * [13] Pedro Larranaga, Borja Calvo, Roberto Santana, Concha Bielza, Josu Galdiano, Inaki Inza, José A Lozano, Rubén Armananzas, Guzmán Santafé, Aritz Pérez, et al.  Machine learning in bioinformatics.  Briefings in bioinformatics, 7(1):86–112, 2006.
  * [14] Belay Enyew Chekol and Hani Hagras.  Employing machine learning techniques for the malaria epidemic prediction in ethiopia.  In 2018 10th Computer Science and Electronic Engineering (CEEC), pages 89–94. IEEE, 2018.
  * [15] Muthoni Masinde.  Africa’s malaria epidemic predictor: Application of machine learning on malaria incidence and climate data.  In Proceedings of the 2020 the 4th International Conference on Compute and Data Analysis, pages 29–37, 2020.
  * [16] Odu Nkiruka, Rajesh Prasad, and Onime Clement.  Prediction of malaria incidence using climate variability and machine learning.  Informatics in Medicine Unlocked, 22:100508, 2021.
  * [17] T Sajana and MR Narasingarao.  An ensemble framework for classification of malaria disease.  ARPN Journal of Engineering and Applied Sciences, 13(9):3299–3307, 2018.
  * [18] Lionel Divin Mfisimana, Emile Nibayisabe, Kingsley Badu, and David Niyukuri.  Exploring predictive frameworks for malaria in burundi.  Infectious Disease Modelling, 7(2):33–44, 2022.
  * [19] Daniel J Stekhoven.  missforest: Nonparametric missing value imputation using random forest.  Astrophysics Source Code Library, 2015.


# Modelling Arbitrary Complex Dielectric Properties – an automated
implementation for gprMax

Sylwia Majchrowska1, Iraklis Giannakis2, Craig Warren3, and Antonios
Giannopoulos4<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>1Department of Optics
and Photonics, Wroclaw University of Science and Technology, Wrocław, Poland
2School of Geosciences, University of Aberdeen, Aberdeen, United Kingdom
3Department of Mechanical and Construction Engineering, Northumbria
University, Newcastle, United Kingdom  4School of Engineering, The University
of Edinburgh, Edinburgh, United Kingdom

###### Abstract

There is a need to accurately simulate materials with complex electromagnetic
properties when modelling Ground Penetrating Radar (GPR), as many objects
encountered with GPR contain water, e.g. soils, curing concrete, and water-
filled pipes. One of widely-used open-source software that simulates
electromagnetic wave propagation is gprMax. It uses Yee’s algorithm to solve
Maxwell’s equations with the Finite-Difference Time-Domain (FDTD) method. A
significant drawback of the FDTD method is the limited ability to model
materials with dispersive properties, currently narrowed to specific set of
relaxation mechanisms, namely multi-Debye, Drude and Lorentz media.
Consequently, modelling any arbitrary complex material should be done by
approximating it as a combination of these functions. This paper describes
work carried out as part of the Google Summer of Code (GSoC) programme 2021 to
develop a new module within gprMax that can be used to simulate complex
dispersive materials using multi-Debye expansions in an automatic manner. The
module is capable of modelling Havriliak-Negami, Cole-Cole, Cole-Davidson,
Jonscher, Complex-Refractive Index Models, and indeed any arbitrary dispersive
material with real and imaginary permittivity specified by the user.

###### Index Terms:

Havriliak-Negami, Cole-Cole, FDTD, Jonsher, multi-Debye, Electrodynamics, GPR

##  I Introduction

gprMax (http://www.gprmax.com) was created in the mid-1990s as a tool that can
be used to model Ground Penetrating Radar (GPR) responses from arbitrarily
complex targets [1]. gprMax was originally written in C programming language,
and became one of the most widely used Finite-Difference Time-Domain (FDTD)
based solvers in the GPR community.

In recent years gprMax has been re-developed and signiﬁcantly modernised to
enable more complex models and advanced features [2]. At present, the open-
source code is a combination of Python and Cython programming languages, which
make it more flexible, user-friendly, and provide the opportunity to
incorporate commonly used Python libraries such as Numpy and SciPy. Moreover,
the utilisation of Cython – a superset of the Python language – allows it to
maintain a similar computation speed to the original C-based code.

The numerical modelling approach employed in gprMax is based on Yee’s [3]
algorithm (with second order accurate derivatives in space and time), which is
used to solve Maxwell’s equations with the FDTD method. However, the FDTD
method is limited regarding modelling dispersive materials. Due to its time-
domain nature, arbitrary dispersive materials cannot be directly implemented
in a straightforward manner [4]. To overcome this, the given dielectric
spectrum is approximated by functions compatible with the time-domain
numerical solver [5]. The most common function employed for this is the multi-
pole Debye expansion. The package introduced in this work incorporates into
gprMax three optimization approaches to fit a multi-pole Debye expansion to
dielectric data.

The paper is organized as follows: Section II provides an overview of the
properties of linear dispersive media, and the relaxation functions that were
implemented; Section III describes the design of the provided package along
with the optimisation methods used; Section IV presents the performance of the
resulting Debye models in terms of computational time and dielectric model
accuracy; and ﬁnally Section V discusses the outcomes of the new implemented
features.

##  II Dielectric Properties of Matter

###  II-A Linear Dispersive Materials

The dielectric properties of linear dispersive materials can be expressed as a
frequency-dependent complex function

$\epsilon(\omega)=\epsilon^{{}^{\prime}}(\omega)-j\epsilon^{{}^{\prime\prime}}(\omega),$
(1) $\mu(\omega)=\mu^{{}^{\prime}}(\omega)-j\mu^{{}^{\prime\prime}}(\omega),$
(2)

where $\omega$ is the angular frequency, $\epsilon^{{}^{\prime}}$ and
$\mu^{{}^{\prime}}$ stand for the real parts of electric permittivity and
magnetic permeability respectively; while $\epsilon^{{}^{\prime\prime}}$ and
$\mu^{{}^{\prime\prime}}$ denote the imaginary parts. In general, as a hard
and fast rule, the real part dictates the velocity of the medium while the
imaginary part is related to the electromagnetic losses [6]. Values of complex
relative permittivity and magnetic permeability can be obtained with a number
of different techniques proposed and developed over the last decades [7, 8,
9].

It appears that the most common group of materials, such as water, soil and
biological tissues, exhibit dispersive behavior over a wide frequency range
and linear behavior for commonly used field intensities [10]. The most
important functions used to describe the behaviour of linear dispersive
materials are the Debye, Lorentz, Drude, Havriliak-Negami [4], Jonsher
equations, and the Complex Refractive Index Model (CRIM) [11].

###  II-B Dielectric Relaxation Mechanisms

Many materials of practical interest in GPR are linear, isotropic, and have
insignificant magnetic response. In such cases they can be modeled using a
scalar complex relative permittivity function. Where frequency dependence of
electric permittivity is required, typically the single-pole Debye equation is
used [6]

$\epsilon(\omega)=\epsilon_{\infty}+\frac{\epsilon_{s}-\epsilon_{\infty}}{1+j\omega\tau_{0}},$
(3)

where $\epsilon_{s}$ and $\epsilon_{\infty}$ are the values of the real part
of the complex relative permittivity at static and infinite frequency,
respectively, and $\tau_{0}$ is the relaxation time, which is directly related
to relaxation frequency – the frequency for which the imaginary part as well
as the first derivative of the real part reach their maximum value (see Fig.
1).

However, the single-pole Debye equation is a poor model of dielectric behavior
for most materials over wide frequency ranges. To overcome this, the
empirically derived Havriliak–Negami relaxation was proposed as a modification
of the single-pole Debye relaxation model. The equation additionally has two
exponential parameters, and is expressed as

$\epsilon(\omega)=\epsilon_{\infty}+\frac{\epsilon_{s}-\epsilon_{\infty}}{\left(1+\left(j\omega\tau_{0}\right)^{\alpha}\right)^{\beta}},$
(4)

where $\alpha$ and $\beta$ are positive real constants ($0\geq\alpha,\beta\leq
1$). From this model, the Cole-Cole equation [12] setting $\beta=1$, and Cole-
Davidson equation [13] for $\alpha=1$ can be derived. The Debye equation is
obtained with $\alpha=1$ and $\beta=1$. The Havriliak-Negami function was
first used to describe the dielectric properties of polymers [14], whereas the
Cole-Cole equation is mainly used to model biological tissues [15] and liquids
[16].

On the other hand, the Jonscher function is mainly used to describe the
dielectric properties of concrete [17] and soils[18]. The frequency domain
expression of the Jonscher function in the constant Q-factor approach (or
quality Q-factor) [19] is given by

$\epsilon(\omega)=\epsilon_{\infty}+A_{p}\left(-j\frac{\omega}{\omega_{p}}\right)^{n_{p}},$
(5)

where $A_{p}$ is a Jonscher parameter (positive real constant), $n_{p}$
characterizes the change in amplitude as a function of frequency (it varies
between 0 for materials with high dielectric loss and 1 for materials without
dielectric losses), and $\omega_{p}$ is the reference frequency, arbitrarily
chosen.

Other types of relaxation functions estimate the bulk permittivity of
heterogeneous materials. Here, CRIM [20] was suggested as an improvement to
the model initially proposed by Brown [21]. The model illustrates the
dielectric properties of the mixture with respect to the dielectric properties
and the volumetric fractions of its components, and it is defined as

$\epsilon(\omega)^{a}=\sum_{i=1}^{m}f_{i}\epsilon_{m,i}(\omega)^{a},$ (6)

where $a$ is a shape factor (commonly set to 0.5), $f_{i}$ and
$\epsilon_{m,i}$ represents volumetric fraction and relative permittivity for
$i_{\mathrm{th}}$ material, respectively. The CRIM has been established as the
mainstream methodology in the GPR community mainly due to its simplicity and
its straightforward implementation.

The Havriliak-Negami and Jonscher functions as well as steady Q-factor, CRIM,
and experimentally obtained data cannot be directly implemented into the FDTD
algorithm. An approach used by many researchers is the approximation of the
complex permittivity using a multi-pole Debye expansion

$\epsilon(\omega)=\epsilon_{\infty}+\sum_{n=1}^{N}\frac{\Delta\epsilon_{n}}{1+j\omega\tau_{0,n}},$
(7)

with $\Delta\epsilon_{n}$ and $\tau_{0,n}$ being the change in permittivity
and relaxation times of the Debye dispersion, respectively. Parameter $N$
indicates the number of Debye poles (the total number of Debye expansion
components).

Figure 1: Both real and imaginary parts of a Debye function with
$\epsilon_{s}=15\epsilon_{0}$ ($\epsilon_{0}$ is the electric permittivity of
free space), $\epsilon_{\infty}=\epsilon_{0}$ and $\tau_{0}=0.1\,\mathrm{ns}$.
The relaxation frequency equals $f_{0}=\frac{1}{2\pi\tau_{0}}$.

##  III Package overview

###  III-A Code structure

The current package features a stochastic global optimisation approach to fit
a multi-pole Debye expansion to dielectric data. The user can choose between
Havriliak-Negami, Jonscher, and CRIM. As mentioned previously, the Havriliak-
Negami relaxation function is an inclusive function that holds as special
cases the widely-used Cole-Cole and Cole-Davidson models. The present package
can be also used to fit arbitrary dielectric data derived experimentally or
calculated using some formulated expression.

The package consists of two main modules:

  * •

Debye_fit – being the main core of the package, containing the definition of
the implemented relaxation functions and methods to run the optimization
procedure,

  * •

optimization – containing the implementation of selected global optimization
approaches.

In the package we used three types of optimisation algorithms for determining
the relaxation times, namely Dual Annealing (DA) [22], Differential Evolution
(DE) [23] and Particle Swarm Optimization (PSO) [24] techniques. To calculate
values of the weights we used the Damped Least-Squares (DLS) method. Our
implementations are mainly based on modules provided in the pyswarm and SciPy
optimize packages, which deliver various methods to minimize objective
functions.

###  III-B Fitting procedure

Accurately fitting Debye coefficients to a formulated relaxation function or a
set of dielectric measurements across a very large frequency-range is
problematic. In literature various fitting methods have been investigated
previously including genetic algorithms [25, 26], and hybrid particle-swarm
least squares [5] optimization procedures. The whole idea is based on
minimization of distance between the calculated values of the relaxation model
and fitted multi-pole Debye expansion at individual frequency points in some
defined range.

The optimisation method employed in the package is a hybrid linear-nonlinear
optimisation approach. As a default we used a slightly adjusted optimization
procedure to approximate the complex permittivity using the multi-pole Debye
function expansion, as proposed by Kelley et. al. [5]. In their work, LS was
used to determine the weights and PSO to determine the relaxation frequencies
($\omega_{0,n}=1/\tau_{0,n}$). Our modification is based on changing the sign
to opposite in case of negative weights in order to introduce a large penalty
in the optimisation process thus indirectly constraining the weights to always
be positive. Furthermore we added the real part to the cost action to avoid
possible instabilities to arbitrary given functions that do not follow the
Kramers–Kronig relationship. These changes overcome some instability issues
and thus make the process more robust and faster.

During creation of a relaxation object the user has a possibility to choose
one from three available optimization techniques (and optionally some
additional fixed parameters needed to completely specify the behavior of the
optimizer) to calculate relaxation times. For the PSO procedure it is possible
to plot the average error between actual and the approximated value of an
objective function during the optimization process (see Fig. 2). Additionally
we give a possibility to automatically set number of Debye poles. The
estimation of the number of Debye poles is done iteratively starting from the
one-pole Debye model. The optimization stops when it reaches the desired
accuracy (error below 5%) or maximum number of iteration (20 Debye poles). The
procedure is launched when the parameter determining the number of Debye poles
is set to -1.

Figure 2: Example of fitting procedure for PSO algorithm.

##  IV Multi-pole Debye models for FDTD simulations

To illustrate the application of the package, a Havriliak-Negami permittivity
model was approximated over the frequency range of 10 MHz to 100 GHz. The
model parameters for the example are $\epsilon_{s}=8.6\epsilon_{0}$,
$\epsilon_{\infty}=2.7\epsilon_{0}$, $\tau_{0}=0.94\,\mathrm{ns}$,
$\alpha=0.91$, and $\beta=0.45$. While the frequency range is wide, the number
of Debye poles was set to 5 as in [5]. The results of the fitting procedure
are presented in Fig. 3.

Figure 3: Real and imaginary parts of the five-pole Debye function expansions
achieved with PSO algorithm compared to those of the Havriliak-Negami
permittivity model [5], and relative approximation error for both real and
imaginary parts.

We compared results from all algorithms in terms of computational speed and
accuracy. The average fractional error of each algorithm has been calculated
separately for both real and imaginary part of the relative permittivity and
then summed. The obtained parameters are presented in Table I. It appears that
the PSO-DLS combination is the fastest optimization procedure, which achieves
relatively good approximation results for the presented example (below 10% of
relative average error).

TABLE I: Summary of results for fitting relative permittivity. Algorithm | Duration (s) | Average Error (%)
---|---|---
PSO-DLS | 0.97 | 8.48
DA-DLS | 2.56 | 4.37
DE-DLS | 3.61 | 4.48

Multi-pole Debye, Drude, and Lorentz functions are already available in gprMax
input files. This made it easier to integrate this package into the software
itself. Listing 1 gives an example of the command to use a described above
example using PSO method inside gprMax software.

Listing 1: A 5-pole Debye expansion to approximate Havriliak-Negami
permittivity model.

⬇

#havriliak_negami: 1e7 1e11 0.91 0.45 2.7 5.9 9.4e-10 0.1 1 0 5 Kelley

The gprMax commands (like $\\#havriliak\\_negami$, $\\#jonscher$, $\\#crim$,
$\\#raw\\_data$) to model dispersive materials deﬁne the basic material
properties along with the relaxation function parameters. In Listing 1 the
parameters for the $\\#havriliak\\_negami$ command are given going from left
to right: the lower and upper frequency bounds (Hz), relaxation function
parameters (here $\alpha$ and $\beta$ as for eq. 4), and basic material
properties as real relative permittivity at infinite frequency, the difference
between the static relative permittivity and the relative permittivity at
infinite frequency, relaxation time (seconds), the conductivity
(Siemens/metre), the relative permeability, the magnetic loss (Ohms/metre),
the number of Debye poles, and an identifier for the material (text label).
Additionally at the end of each command user can put optional integer, which
controls the seeding of the random number generator used in stochastic global
optimizer. During fitting the hybrid PSO algorithm is used as a default. Full
description about other hashtag commands could be found in the documentation.

Listing 2 gives an example of the output for the $\\#havriliak\\_negami$
command described above. Line 1 defines the basic material properties with the
$\\#material$ command, and in line 2 the $\\#add\\_dispersion\\_debye$ command
adds dispersive behaviour to the material based on the Debye formulation. In
the definition of the $\\#add\\_dispersion\\_debye$ command the number of
Debye poles, and then the the difference between the static relative
permittivity and the relative permittivity at infinite frequency for the
$n_{\mathrm{th}}$ Debye pole, the relaxation time (seconds) for the
$n_{\mathrm{th}}$ Debye pole are given. In case of running stand alone script
these commands supposed to be added to the input file. On the other hand, the
use of the $\\#havriliak\\_negami$ (Listing 1) command will pass the
approximated parameters inside gprMax automatically.

Listing 2: Example output for 5-pole Debye expansion to approximate Havriliak-
Negami permittivity model.

⬇

#material: 2.8345 0.1 1 0 Kelley

#add_dispersion_debye: 5 2.3563 4.3677e-10 0.6736 1.1623e-11 1.1009 1.5180e-9
1.2926 9.1048e-11 0.3091 1.1131e-12 Kelley

##  V Conclusions

The current work describes the addition of a new advanced modelling feature in
gprMax capable of simulating complex dielectric properties defined by some
common functions such as Havriliak-Negami, Jonscher, CRIM, and any arbitrary
user-specified dispersive dielectric properties. The hybrid linear-nonlinear
optimisation approach employed here provides an effective, and accurate
optimisation procedure to fit a multi-pole Debye expansion to the given
dielectric data and subsequently implement it to gprMax in an automatic
manner.

## Acknowledgment

The project was funded via Google Summer of Code (GSoC) 2021 programme. GSoC
initiative is a global program focused on bringing more student developers
into open source software development.

## References

  * [1] A. Giannopoulos, “Modelling ground penetrating radar by gprmax,” Construction and building materials, vol. 19, no. 10, pp. 755-–762, August 2005.
  * [2] C. Warren, A. Giannopoulos, and I. Giannakis, “gprMax: Open source software to simulate electromagnetic wave propagation for Ground Penetrating Radar,” Computer Physics Communications, vol. 209, pp. 163–170, September 2016.
  * [3] K. S. Yee, “Numerical solution of initial bound-ary value problems involving maxwells equations in isotropic media, Antennas and Propagation,” IEEE Transactions, vol. 14, no. 3, pp. 302-–307, May 1966.
  * [4] I. Giannakis and A. Giannopoulos, “A Novel Piecewise Linear Recursive Convolution Approach for Dispersive Media Using the Finite-Difference Time-Domain Method,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 5, pp. 2669-2678, May 2014
  * [5] D. F. Kelley, T. J. Destan and R. J. Luebbers, “Debye Function Expansions of Complex Permittivity Using a Hybrid Particle Swarm-Least Squares Optimization Approach,” IEEE Transactions on Antennas and Propagation, vol. 55, no. 7, pp. 1999–2005, July 2007.
  * [6] C. A. Balanis, “Advanced Engineering Electromagnetics,” New York:Wiley, Sec. 2.8.1, 1989
  * [7] M. Afsar, J. Birch, R. Clarke, and G. Chantry, “The measurement of the properties of materials,” Proceedings of the IEEE, vol. 74, no. 1, pp. 183-–99, 1986.
  * [8] P. Queffelec, P. Gelin, J. Gieraltowski, and J. Loaec, “Microstrip device for the broad band simultaneous measurement of complex permeability and permittivity,” IEEE Transactions on Magnetics, vol. 30, no. 2 pt 1, pp. 224–231, 1994.
  * [9] J. Baker-Jarvis, C. Jones, B. Riddle, M. Janezic, R. Geyer, Jr. J. H. Grosvenor, and C. Weil, “Dielectric and magnetic measurements: a survey of nondestructive, quasi-nondestructive, and process-control techniques,” Research in Nondestructive Evaluation, vol. 7, no. 2–3, pp. 117-–36, 1995.
  * [10] A. Taflove, and S. Hagness, “Computational Electromagnetics: The Finite-Difference Time-Domain Method,” Artech House Norwood, MA, 3rd edition, 2005.
  * [11] H. Zadhoush, A. Giannopoulos, I. Giannakis, “Optimising the Complex Refractive Index Model for Estimating the Permittivity of Heterogeneous Concrete Models,” Remote Sens. vol. 13, no. 723, 2021.
  * [12] K. S. Cole, and R. H. Cole, “Dispersion and absorption in dielectrics, alternating current characteristic,” Journal of Chemical Physics, vol. 9, pp. 341–351, 1941.
  * [13] M. F. Manning and M. E. Bell, “Electrical conduction and related phenomena in solid dielectrics,” Rev. Mod. Phys., vol. 12, pp. 215–257, 1940.
  * [14] S. Havriliak, and S. Negami, “A complex plane representation of dielectric and mechanical relaxation processes in some polymers,” Polymer, vol. 8, pp. 161–210, 1967.
  * [15] D. Ireland, A. Abbosh, “Modeling Human Head at Microwave Frequencies Using Optimized Debye Models and FDTD Method,” IEEE Trans. Antennas and Propagation, vol. 61, no. 4, pp. 2352–2355, 2013.
  * [16] J. Barthel, R. Buchner, and M. Münsterer, “Electrolyte Data Collection, Part 2: Dielectric Properties of Water and Aqueous Electrolyte Solutions,” ser. Chemistry Data Series. Frankfurt am Mein, Germany:DECHEMA, vol. XII, 1995.
  * [17] T. Bourdi, J. Rhazi, F. Boone, and G. Ballivy, “Application of Jonscher model for the characterization of the dielectric permittivity of concrete,” Journal of Physics D: Applied Physics, vol. 41, pp. 205410, October 2008.
  * [18] E. Kjartansson, “Constant Q-wave propagation and attenuation,” Journal of Geophysical Research, vol. 84, pp. 4737–4748, 1979.
  * [19] M. Bano, “Modelling of GPR waves for lossy media obeying a couplex power law of frequency for dielectric permittivity,” Geophysical Prospecting, vol. 41, pp 11–26, 2004.
  * [20] J.R. Birchak, C.G. Gardner, J.E. Hipp, J.M Victor, “High dielectric constant microwave probes for sensing soil moisture,” Proc. IEEE, vol. 62, pp. 93–-98, 1974.
  * [21] W. F. Brown, “Dielectrics,” Encyclopedia of Physics, vol. 17, Berlin: Springer, 1956.
  * [22] C. Tsallis, DA Stariolo, “Generalized Simulated Annealing,” Physica A, vol. 233, pp. 395–406, 1996.
  * [23] R. Storn, and K. Price, “Differential Evolution - a Simple and Efficient Heuristic for Global Optimization over Continuous Spaces,” Journal of Global Optimization, vol. 11, 341–359, 1997.
  * [24] J. Kennedy and R. Eberhart, “Particle swarm optimization,” Proceedings of ICNN’95 - International Conference on Neural Networks,” vol.4, pp. 1942–1948, 1995.
  * [25] J. Clegg, M. P. Robinson, “A genetic algorithm for optimizing multi-pole Debye models of tissue dielectric properties,” Physics in Medicine and Biology, vol. 57, pp. 6227-–6243, 2012.
  * [26] F. Krewer, Morgan, M. O’Halloran, “Development of accurate multi-pole Debye functions F. for electromagnetic tissue modelling using a genetic algorithm,” Progress In Electromagnetics Research Letters, vol. 43, pp. 137–-147, 2013.