Yu and Haskins

# KNN, An Underestimated Model for Regional Rainfall Forecasting

Ning Yu
Department of Computing Sciences
SUNY Brockport
350 New Campus Drive
Brockport  NY  USA 14420
<EMAIL_ADDRESS>

Timothy Haskins
Department of Computing Sciences
SUNY Brockport
350 New Campus Drive
Brockport  NY  USA 14420
<EMAIL_ADDRESS>Roberto Szechtman
Han L. Yooh
Jonathan Harker
Department of Operations Research
Naval Postgraduate School
1 University Circle
Monterey  CA  USA
{scs21hanlee  scs21szechtman
<EMAIL_ADDRESS>

Enlu Zhou

H. Milton School of Industrial & Systems Engineering
Georgia Institute of Technology
755 Ferst Drive NW
Atlanta  GA  USA
<EMAIL_ADDRESS>

## Abstract

Regional rainfall forecasting is an important issue in hydrology and
meteorology. This paper aims to design an integrated tool by applying various
machine learning algorithms, especially the state-of-the-art deep learning
algorithms including Deep Neural Network, Wide Neural Network, Deep and Wide
Neural Network, Reservoir Computing, Long Short Term Memory, Support Vector
Machine, K-Nearest Neighbor for forecasting regional precipitations over
different catchments in Upstate New York. Through the experimental results and
the comparison among machine learning models including classification and
regression, we find that KNN is an outstanding model over other models to
handle the uncertainty in the precipitation data. The data normalization
methods such as ZScore and MinMax are also evaluated and discussed.

Keywords: rainfall forecasting, k-nearest neighbor, deep and wide neural
network, reservoir computing, long short term memory.

##  1 Introduction

New York historically had sufficient precipitation until recently, with
intense drought occurring over the 2016 growing season, especially in western
New York [Todaro (2018)]. The observed precipitation in 2016 was less than
normal, with shortfalls of 4-8 inches being common in the 90 days leading up
to the drought watch. Accurate rainfall forecasting is important for planning
in agriculture and other relevant activities.

Although a number of modern algorithms and applications have been used to
forecast rainfall, there are two categories of approaches to solve the
problem. The first category models the underlying physical principles for the
rainfall process. However, it is thought not feasible limited by the complex
climatic system in various spatial and temporal dimensions. A second category
is based on the data mining and pattern recognition, which attempts to mine
rainfall patterns and learn the knowledge from numerous features and a large
volume of data. Historical meteorological data including precipitation data
are used to feed and train the recognition model and further predict the
evolution of other storms. In the recent years, with the advance of big data
and deep learning technology, this type of approach has drawn more attentions
from researchers. It does not require a thorough understanding of the climatic
laws but it does need the data modeling for data mining and pattern
recognition [Luk et al. (2001)].

There are a number of reported research that have used machine learning
algorithms to solve problems in hydrology. These modern machine learning
algorithms include Deep Neural Network (DNN), Wide Neural Network (WNN), Deep
and Wide Neural Network (DWNN), Reservoir Computing (RC), Long Short Term
Memory (LSTM), Support Vector Machine (SVM), K-Nearest Neighbor (KNN), and so
forth.

As a type of Artificial Neural Networks (ANN), DNN can emulate the process of
the human nervous system and have proven to be very powerful in dealing with
complicated problems, such as computer vision and pattern recognition.
Moreover, DNN is computationally robust even when input data contains lots of
uncertainty such as errors and incomplete, as such examples are very common in
rainfall prediction.

A pioneer two-dimensional ANN model [French et al. (1992)] has been considered
as the first Neural Network (NN) application to simulate complex geophysical
processes and forecast the rainfall, although it was limited by many aspects
including insufficient neural network configurations and a mathematical
rainfall simulation model that was used to generate the input data. More
applications of NN model were developed later. [Koizumi (1999)] utilized an
ANN model for radar, satellite and weather-station data together to provide
better performance than the persistence forecast, the linear regression
forecasts, and the numerical model precipitation prediction. By comparing
several short-time rainfall prediction models including linear stochastic
auto-regressive moving-average (ARMA) models, artificial neural networks (ANN)
and the non-parametric nearest-neighbours method, [Toth et al. (2000)]
revealed that those ANN models provide a significant improvement in real-time
flood forecasting. [Luk et al. (2001)] developed three types of Artificial
Neural Networks including time delay neural networks for rainfall forecasting.
They revealed that the rainfall time series have very short term memory
characteristics.

[Hung et al. (2009)] applied ANN model for real time rainfall forecasting and
flood management in Bangkok, Thailand. The developed ANN model was fed with
meteorological parameters (relative humidity, air pressure, wet bulb
temperature and cloudiness), the rainfall at the point of forecasting and
rainfall at the surrounding stations that fueled the ANN model to predict the
rainfall at any moment.

In recent years, the performance of three Artificial Neural Networks (ANN)
approaches were evaluated in the forecasting of the monthly rainfall anomalies
for Southwestern Colombia [Canchala et al. (2020)]. A forecasting performance
was found suitable for the regional rainfall forecasting in Southwestern
Colombia and the combination of ANN approaches had demonstrated the ability to
provide useful prediction information for the decision-making.

Google researchers have invented Wide & Deep Neural Network (WDNN) model for
recommender systems originally, where it contains two components: wide neural
network and deep neural network. The wide neural network (WNN) was designed to
effectively memorize the feature interactions through a wide set of cross-
product feature transformations. The deep neural networks (DNN) were used to
generalize better to unseen feature combinations through low-dimensional dense
embeddings learned for the sparse features [Cheng et al. (2016)]. WDNN has
been used in various applications, such as clinical regression prediction for
Alzheimer patients [Pölsterl et al. (2019)]. In the latest publication,
[Bajpai and Bansal (2021)] analyzed and evaluated three deep learning
approaches, one dimensional Convolutional Neutral Network, Multi-layer
Perceptron and Wide Deep Neural Networks for the prediction of summer monsoon
rainfall in Indian state of Rajasthan and found that Deep and Wide Neural
Network (DWNN) can achieve the best performance over the rest two deep
learning methods.

As a type of recurrent neural network (RNN), Long Short-Term Memory neural
network (LSTM) has been studied by many researchers in recent literature. [Liu
et al. (2020)] used LSTM to simulate rainfall–runoff relationships for
catchments with different climate conditions. The LSTM model then was coupled
with the k-nearest neighbor (KNN) algorithm as an updating method to further
boost the performance of LSTM method. LSTM-KNN model was validated by
comparing the single LSTM and the recurrent neural network (RNN). [Chhetri et
al. (2020)] studied LSTM, Bidirectional LSTM (BLSTM), and Gated Recurrent Unit
(GRU), another variant of RNN. LSTM recorded the best Mean Square Error (MSE)
score outperformed other six different existing models including Linear
Regression, Multi-Layer Perceptron (MLP), Convolutional Neural Network (CNN),
GRU, and BLSTM. They further modified the model by combining LSTM and GRU to
carry out monthly rainfall prediction over a region in the capital of Bhutan.

As one of machine learning models, K-nearest neighbor model has shown a
promising performance in climate prediction. [Jan et al. (2008)] applied KNN
for climate prediction by using the historical weather data of a region such
as rain, wind speed, dew point, temperature, etc. [Huang et al. (2017)]
developed a K-nearest neighbor (KNN) algorithm to offer robustness the
irregular class distribution of the precipitation dataset and made a sound
performance in precipitation forecast. However, they did not show the
comparison with other advanced machine learning algorithms such as deep
learning, limiting its application on the rainfall forecasting. [Yang et al.
(2020)] developed an Ensemble–KNN forecasting method based on historical
samples to avoid uncertainties caused by modeling inaccuracies. [Hu et al.
(2013)] proposed a model that combined empirical mode decomposition (EMD) and
the K-nearest neighbor model to improve the performance of forecasting annual
average rainfall. [Wu (2009)] proposed a KNN nonparametric regression model
used to improve an ANN-based model.

Like multilayer perceptrons and radial basis function networks, support vector
machines can be used for pattern classification and nonlinear regression. SVM
is found to be a significant technique to solve many classifications problem
in the last couple of years. [Hasan et al. (2015)] exhibited a robust rainfall
prediction technique using Support Vector Regression (SVR). Support Vector
Machine based approaches also have illustrated the quite capability for
rainfall forecasting. Support vector regression combined with Singular
Spectrum Analysis (SSA) has shown its efficiency for monthly rainfall
forecasting in two reservoir watersheds [Bojang et al. (2020)].

The Reservoir Computing framework has been known for over a decade as a state-
of-the-art paradigm for the design of RNN [Coble and Yu (2020)]. Among the
models of RC instances, the echo state network (ESN) represents a type of the
most widely known schemes with a strong theoretical support and various
applications [Gallicchio and Micheli (2017)]. [Yen et al. (2019)] used the ESN
and the DeepESN algorithms to analyze the meteorological hourly data in Taiwan
and showed that DeepESN was better than that by using the ESN and other
neuronal network algorithms. Backpropagation Algorithms and Reservoir
Computing in Recurrent Neural Networks were adopted to predict the complex
spatiotemporal dynamics [Vlachas et al. (2020)].

The precipitation prediction models in this study are based on not only the
implementation of the state-of-the-art machine learning algorithms but also
the design of an integrated tool for these machine learning algorithms to
provide better regional rainfall forecasting. A variety of machine learning
algorithms are applied and compared including DNN, WNN, DWNN, RC, LSTM, SVM,
KNN for forecasting regional precipitations over different catchments in
Upstate New York. These algorithms were applied to forecast rainfall at
multiple locations simultaneously for the next hour based on previously
observed precipitation data. The important issues of data normalization and
comparison of algorithms are discussed also.

##  2 Models and Algorithms

A set of algorithms of DNN, WNN, DWNN, RC, LSTM, SVM, KNN are used to forecast
regional rainfall. Both classification and regression models of these
algorithms are implemented in our integrated tool. For the classification
models, the threshold of rainfall and non-rainfall is 0.01 inch. Any of below
the threshold is defined as very light precipitation according to the
categories in hydrology. For the regression models, the root of mean-squared
error (RMSE) and mean-squared error (MSE) are used to measure the performance
of these regression models. We describe these main models and algorithms used
in our study as following subsections.

###  2.1 Wide Neural Network

The core of a wide neural network is a generalized linear model as shown in
(1).

$y=w^{T}x+b,$ (1)

where $x$ is a vector of $d$ features, $w$ is the weights or parameters, and
$b$ is the bias. The $d$ feature set may be raw input features or transformed
features. A cross-product transformation is utilized to capture the
interactions between the binary features, which adds a non-linear factor to
the generalized linear model, as shown in (2).

$\phi_{k}(x)=\prod_{i=1}^{d}x_{i}^{c_{ki}},c_{ki}\in\\{0,1\\},$ (2)

where $c_{ki}$ is 1 if the $i$-th feature is part of the $k$-th transformation
$\phi_{k}$ and 0 if the $i$-th feature is not part of the $k$-th
transformation $\phi_{k}$. The generalized linear models with cross-product
feature transformations can memorize these exception rules with much fewer
parameters [Cheng et al. (2016)].

###  2.2 Deep Neural Network

Figure 1: Structure of perceptron

Artificial neural Networks are generalized models of biological neuron system.
Deep neural network is a type of artificial neural networks, which contains a
number of layers and can massively distribute processing over layers of
neurons or perceptrons. Each neuron or perceptron is the elementary unit
containing a simple linear summing function and an activation function as
shown in Figure 1. The summing function is illustration in Equation 3:

$s_{k}=\sum_{j=1}^{m}w_{kj}x_{j}+b_{k},$ (3)

where $k$ denotes the identifier of the $k$-th layer, $x_{j}$ is the input of
node $j$, and $b_{k}$ is the bias of layer $k$. The sigmoid function is often
used as an activation function as shown in Equation 4:

$y_{k}(s)=\varphi(s)=\frac{1}{1+e^{-s}}.$ (4)

###  2.3 Wide and Deep Neural Network

Figure 2: An illustration of DWNN

The Wide and Deep Neural Network is composed of a wide neural network and a
deep neural network using a weighted sum of their output log odds as the
prediction and a common logistic loss function for a joint training. In the
joint training, the wide component uses a small number of cross-product
feature transformations to complements the weaknesses of the deep neural
network [Cheng et al. (2016)].

The deep component is a feed-forward neural network, as shown on the right
side of Figure 2 while the wide component is a linear model based neural
network, as shown on the left side of Figure 2. The low-dimensional vectors
are fed into the hidden layers of a neural network in the forward pass
according to Equation 5:

$y_{k+1}=\varphi(w_{k}y_{k}+b_{k}),$ (5)

where $k$ is the layer number and $\varphi$ is the activation function,
$y_{k}$, $w_{k}$, and $b_{k}$ are the output, model weights, and bias at the
$k$-th layer.

###  2.4 K-Nearest Neighbors

KNN algorithm is a non-parametric method that can approximate its outcome by
averaging the observations in the same neighbourhood. The neighbourhood is
determined by the distance between the observations and the input point. The
distance can be measured by many methods such as Manhattan distance, Euclidean
distance, and Minkowski distance. Euclidean distance is calculated as the
square root of the sum of the squared differences between a new point (x) and
an existing point (y) to calculate the distance between the new point and each
training point. Equation 6 shows the distance measurement:

$D(x,y)=\sqrt{\sum_{i=1}^{n}(x_{i}-y_{i})^{2}},$ (6)

where $i$ is the $i$-th feature of nodes $x$ and $y$. KNN can be used for both
classification and regression. In kNN regression, the output value is the
average of the values of k nearest neighbors. In kNN classification, the
output value is the most common class among k nearest neighbors. The value of
k needs to be set manually or can be chosen or using cross-validation to
select the size that minimizes the mean-squared error or maximize the
accuracy.

###  2.5 Reservoir Computing

Figure 3: An illustration of reservoir computing

The reservoir computing (RC) framework consists of three layers: an input
layer, reservoir layer, and output layer as depicted in Figure 3.

The input layer consists of the input weights to each reservoir perceptron. It
is a randomly generated matrix of size ($\\#$ of properties$)\times(\\#$ of
perceptrons), with values uniformly distributed between 0 and $\delta$, where
$\delta$ is a scaling factor. This scaling factor allows us to adjust the
reservoir computing method to training data of various sizes. The input layer
matrix is denoted as $W_{in}$.

The reservoir layer is where the actual RNN is held. The current, $i$-th,
state of the reservoir is described as a length ($\\#$ of perceptrons) vector,
referenced as $\vec{r}(i)$. The perceptron weights are described by a
connectivity/adjacency matrix, $A$. This matrix is an Erdös-Rényi randomly
generated matrix. The eigenvalues of the matrix are normalized, and scaled by
a factor, $\rho$, known as the spectral radius. The spectral radius of the
connectivity matrix determines how the reservoir adapts to changing dynamics.
Usually, a spectral radius of around 0.9 is used in dynamics-prediction
applications. The activation function for the reservoir is given in Equation
7.

$\vec{r}(i)=\tanh[\vec{r}(i-1)\times A+\vec{u}(i)\times W{{}_{in}}],$ (7)

where $\vec{r}(i-1)\times A$ represents the reservoir feedback and is what
allows the reservoir to adapt to new input vectors. The $\vec{u}(i)\times
W_{in}$ term represents the input to the reservoir. The input term and
activation function used varies greatly between different implementations of
reservoir computing [Coble and Yu (2020)].

The output layer is how a prediction is obtained from the reservoir. The
output layer is a ($\\#$ of perceptrons) $\times$( $\\#$ of classes) matrix,
referred to as $W_{out}$. The output is obtained by Equation 8.

$output=\vec{r}(i)\times W{{}_{out}}$ (8)

###  2.6 Support Vector Machine

No matter it is for classification or regression, support vector machine (SVM)
or support vector regressor (SVR) aims to find out the hyperplane. The
difference between SVM and SVR is that one is to segregate the nodes for
classification while another is to find the decision boundary along with
hyperplane and have the least error rate to fit the model. The hyperplane in a
SVM is shown in Equation 9.

$y=\vec{w}\cdot x+b$ (9)

The learning of the hyperplane is done by training the samples using some
kernel. Equation 10 shows the prediction of a new input using a kernel.

$f(x)=\sum_{i}a_{i}K(x_{i},x)+b_{0},$ (10)

where $K$ is the kernel function, $x$ is the input, $x_{i}$ the support
vector, the coefficients $b_{0}$ and $a_{i}$ are estimated from the training
data through the learning algorithm. A common kernel function $K$ can be
linear such as dot product, or non-linear such as Radial basis function kernel
(rbf) as shown in Equation 11,

$K(x_{i},x)=e^{-\frac{\left\|x-x_{i}\right\|^{2}}{2\delta^{2}}},$ (11)

where $\delta$ is a free parameter, and $\left\|x-x_{i}\right\|^{2}$ is the
Euclidean distance.

###  2.7 Long Short Term Memory

Figure 4: An illustration of LSTM structure

The LSTM network accepts the output of the previous moment, the current system
state, and the current system input, updates the system status through the
input gate, the forget gate, and the output gate, then finally outputs the
result [Hochreiter and Schmidhuber (1997)]. As shown in Figure 4, the forget
gate is $f(t)$ , the input gate is $i(t)$, the output gate is $o(t)$, and the
cell state is $C_{t}$.

The forget gate decides the information forgotten and controls the transfer of
information from the previous time as shown in Equation 12:

$f_{t}=\delta(U_{f}x_{t}+W_{f}h_{t-1}).$ (12)

The input gate determines the current time system input, which acts in front
of the forget gate to supplement the latest memory of the network, as shown in
Equation 13:

$i_{t}=\delta(U_{i}x_{t}+Wih_{t-1}).$ (13)

The candidate layer $\overline{C}$ is shown in Equation 14:

$\overline{C_{t}}=tanh(U_{c}x_{t}+W_{c}h_{t-1}).$ (14)

The current cell memory $C_{t}$ is defined in Equation 15:

$C_{t}=f_{t}C_{t-1}+i_{t}\overline{C_{t}}.$ (15)

The output gate $O_{t}$ is shown in Equation 16:

$O_{t}=\delta(U_{o}x_{t}+W_{o}h_{t-1}).$ (16)

The current cell output $h_{t}$ is shown in Equation 17:

$h_{t}=tanh(C_{t})O_{t}.$ (17)

Bidirectional LSTM is viewed as a variant of LSTM. Using bidirectional can
receive the input from two directions: one is from past to future and one is
from future to past that is different from unidirectional LSTM. Gated
Recurrent Unit (GRU) is similar to a LSTM model but it lacks a output gate so
that it has fewer parameters than LSTM. It has demonstrated better performance
in smaller dataset [Cho et al. (2014)].

##  3 Data Preparation

The data was obtained from Iowa State University’s Iowa Environmental Mesonet
which has a portal to download Automated Surface Observing System (ASOS) data
from weather stations around the US. Four regional data sets are collected.
Rochester data came from the ROC weather station. Buffalo data came from the
BUF and IAG weather stations. Syracuse data came from the SYR and FZY weather
stations. Albany data came from the ALB and SCH weather stations. With this
tool, we were able to collect 11 years of historical hourly weather data from
the locations we needed for our project. The features of the data that we
collected for each location were; Hourly precipitation in inches, Temperature
in Fahrenheit, Dew Point in Fahrenheit, Relative Humidity as a percent, Wind
Direction in degrees from North, Wind Speed in knots, Pressure altimeter in
inches, Sea Level Pressure in millibar, and Visibility in miles. The target
was the Precipitation in inches for the next hour in Rochester.

The plot matrix shown in Figure 5 shows the relationships between the numeric
features for Rochester. There are 2 sets of features with clear linear
correlations, that being Temperature and Dew Point as well as Pressure and
Pressure at Sea Level. Temperature and Dew point are commonly used in
calculating Heat Index so it makes sense that they are strongly related. The
two different measurements of air pressure seem nearly identical but as you
can see in the chart, there are outliers showing that they are not exactly the
same. It is for this reason that we decided to keep both of them. Besides the
two sets of strongly correlated features, the other features present a more
complicated relationship. In this plot matrix, the data points in which the
following day had precipitation above 0 inches are marked in green while the
rest of the data points are blue.

Figure 5: The feature matrix in ROC dataset

The raw data was grouped into hour intervals starting with 0 on January 1,
2010, at midnight. These groups were averaged into a single row for all of the
features except p01i which is specifically the accumulated precipitation
within an hour so the max was taken for each group for the p01i feature. All
missing p01i values were set to 0. The decision was made after looking at
random hour intervals whose p01i value was missing and cross-referencing that
time frame with Weather.com. In all instances, Weather.com had 0 inches of
accumulated precipitation where the ASOS p01i values were missing. For the
rest of the features, missing values were interpolated from their previous and
future hour interval values.

From this cleaning and grouping, two separate datasets were created. One
dataset had only data from Rochester. The second dataset had data from
Rochester as well as corresponding data from Buffalo, Syracuse, and Albany.
The second dataset was created by performing a left join of the Rochester data
with a union of the rest of the data. Intervals where Rochester had data but
the others did not, had their values interpolated by the same method that was
used to previously remove missing values.

At this point in the data’s life, two separate versions of the two datasets
were made. One for classification where an interval’s target was a binary 0/1
of whether the next interval’s p01i value exceeded 0.01 inches. Another for
the regressors where an interval’s target was the next interval’s p01i value.
The threshold of 0.01 inches was chosen because although hourly rainfall in
inches is already a very small amount, eliminating the smallest of the values
reduced noise.

For each of the 4 datasets, a min-max normalization and a z-score
normalization version were created. This resulted in 6 datasets for
classification and 6 datasets for regressors. For all of the classifiers and
regressors, the 11 years of hourly data was split into 70% for training and
30% for testing. For all except LSTM, the split was made by shuffling the data
by a given random state and split after shuffling. In the cases where random
states were used to shuffle the data, the random states considered were None,
0, and 42.

##  4 Experimental Results

###  4.1 Classification Models

####  4.1.1 KNN

Each of the 6 datasets was crossed with the 3 different random seeds for a
total of 18 forms of data. Each of those 18 forms was iterated through a
SKLearn KNeighborClassifier 42 times to find the best n_neighbor value. After
756 iterations the best accuracy was 96.34% which was obtained from the
Rochester, Buffalo, Syracuse, and Albany dataset, normalized by z-score, a
random state of None, and 11 n_neighbors.

####  4.1.2 RCC

Reservoir Computing Classifier was investigated because, in a previous
investigation of daily values, RCC had the highest accuracy. For RCC, each of
the 6 datasets was crossed with the same 3 random states as before and crossed
again with 6 different sizes of the reservoir; 50, 100, 200, 400, 600, 1000.
The result of these crosses was 108 unique data forms to train with. The best
accuracy was 95.88% which was obtained from the Rochester only dataset,
normalized by z-score, a random state of None, and a reservoir size of 1000.

####  4.1.3 SVM

There were 72 unique Support Vector Machine Classifiers by training against
the cross of 6 premade datasets, 4 Kernels, and the same 3 Random States as
before. The best accuracy for SVM, and ultimately the best accuracy for all
classifiers, was 96.39% obtained from the ROC+BUF+SYR+ALB, normalized with
ZScore, using None as Random State and the rbf kernel.

####  4.1.4 DNN

72 DNN Classifiers were tested and used by training against the cross of the 6
premade datasets, 4 different amounts of hidden layers, and the 3 Random
States used from before. The 4 different hidden layers were: 2 layers, 3
layers, 4 layers, and 5 layers. Each of the hidden layers had the number of
inputs equal to the number of features which for the ROC datasets was 9 and
for the ROC+BUF+SYR+ALB there were 36 inputs at each hidden layer. The best
accuracy for DNN was 94.88% obtained from the ROC+BUF+SYR+ALB, normalized with
ZScore, using None as Random State and with 5 hidden layers.

####  4.1.5 WNN

There were 18 unique WNN Classifiers by training against the cross of the 6
premade datasets and the 3 Random States used from before. The best accuracy
for WNN was 96.06% obtained from the ROC+BUF+SYR+ALB, normalized with ZScore,
using None as Random State.

####  4.1.6 DWNN

There were 72 unique DWNN Classifiers by training against the cross of the 6
premade datasets, 4 different amounts of hidden layers for the deep aspect,
and the 3 Random States used from before. The 4 different hidden layers were;
2 layers, 3 layers, 4 layers, and 5 layers. Each of the hidden layers had the
number of inputs equal to the number of features which for the ROC datasets
was 9 and for the ROC+BUF+SYR+ALB there were 36 inputs at each hidden layer.
The best accuracy for DWNN was 95.91% obtained from the ROC+BUF+SYR+ALB,
normalized with ZScore, using None as Random State and with 4 hidden layers.

####  4.1.7 LSTM

There were 12 unique LSTM Classifiers by training against the cross of the 6
premade datasets, and 2 different sequence lengths. Random states were not
used for LSTM as it depends on chronologically consecutive rows of data. The
two sequence lengths considered were 3 previous rows of data and 7 previous
rows of data. The best accuracy for LSTM was 94.82% obtained from the
ROC+BUF+SYR+ALB, normalized with MinMax, using a sequence length of 3.

####  4.1.8 Overall Classification Results

Overall, both SVM and KNN classifiers performed the best. Figure 6 (A) and (D)
shows all of the best accuracies for each of the classification methods used
in the mixed and ROC datasets. The plotting for this ranking was zoomed into
the accuracy range of 94.5%-96.5% to better show the differences between the
methods. A further examination of the of the ranking with the normalization
data is shown in Figure 6 (C) and (F). In almost every case except for LSTM,
the best normalization method was ZScore. It can also be seen that in many
cases the non-normalized version of the data performed measurably worse than
its ZScore counterpart. Figure 6 (B) and (E) shows another way to see that
same data that shows more clearly the performance in regards to normalization.

Figure 6: (A)-(C) Classification ranking in the Mixed dataset. (D)-(F)
Classification ranking in the ROC dataset. (A) and (D) The best performance
ranking mixing the normalization. (B) and (E) Accuracy with normalization. (C)
and (F) An overall ranking among models and normalizations.

###  4.2 Regression Models

Being that accuracy is not measured the same for regressors the metrics
recorded were $R^{2}$, Mean Square Error (MSE), Root Mean Square Error (RMSE),
and Pearson Correlation Coefficient (PCC).

Some regressors used Numpy’s Random which was creating results that were not
reproducible. By setting a Numpy Random Seed prior to training, results became
reproducible. A quick and simple investigation was done by iterating through
1-100 for Numpy’s Random Seed. Each iterations seed was used to train and
evaluate a simple Linear Regression model. The results of that investigation
were that a Numpy Random Seed of 46 resulted in the best $R^{2}$ value for a
Linear Regression Model. Note that this is separate from the random states
used in training methods.

The data used to train and test the regressors in this section was a subset of
the data used for the classifiers such that all rows targets were above 0.01,
that target being the next hours p01i value (Precipitation inches). This
decision was made after an initial investigation with the full dataset
resulted in much poorer performances. This was likely due to the overwhelming
ration of 0 precipitation hours compared to non-0 precipitation hours.

####  4.2.1 KNN Regressor

Each of the 6 datasets was crossed with the same 3 random states as before to
make 18 unique forms of data and then each of those 18 was run through a KNN
Regressor 1000 times to find the n_neighbor parameter that resulted in the
best $R^{2}$ metrics. The dataset that had the best $R^{2}$, 0.2095615, was
the Rochester, Buffalo, Syracuse, and Albany dataset, normalized using ZScore,
using None for Random State, and with an n_neighbors value of 42. This dataset
also had the best Pearson Correlation score, 0.46443163. The dataset that had
the best MSE and RMSE, 0.005347 and 0.073122, was the Rochester dataset,
normalized using MinMax, using Random State of None, and a n_neighbor value of
51. Table 1 shows the top 2 best metric performers in KNN regressor.

Table 1: Top 2 best metric performers in KNN regressor. Dataset | Normalization | Random | Neighbor | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---|---
Mixed | ZScore | none | 42 | 0.209561502 | 0.006806884 | 0.082503842 | 0.464431631
Mixed | MinMax | None | 38 | 0.164908114 | 0.00707709 | 0.084125443 | 0.418229929
ROC | MinMax | None | 51 | 0.168353653 | 0.005346776 | 0.07312165 | 0.411062473
ROC | None | None | 65 | 0.106086597 | 0.006163514 | 0.078508053 | 0.326574449

####  4.2.2 Linear Regressor

Each of the 6 datasets were crossed with the same 3 random states as before
and run through a linear regression model. The best metrics were obtained from
the Rochester dataset, not normalized, and a random state of None. Those
values were $R^{2}$ = 0.178324, Pearson Correlation = 0.422666, MSE =
0.006913, and RMSE = 0.083146. Table 2 shows the top two performing Linear
Regression Models and their respective performance metrics.

Table 2: Top 2 best metric performers in KNN regressor. Dataset | Normalization | Random | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---
Mixed | ZScore | none | 0.1464424 | 0.007025658 | 0.083819199 | 0.388019644
Mixed | MinMax | None | 0.13850212 | 0.008926641 | 0.094480901 | 0.375745596
ROC | None | None | 0.178324011 | 0.006913207 | 0.083145696 | 0.422665921
ROC | ZScore | None | 0.143381358 | 0.009098939 | 0.095388361 | 0.380159271

####  4.2.3 SVM Regressor

Each of the 6 datasets was crossed with the same 3 random states as before and
crossed again with the 4 different kernel options: linear, poly, rbf, and
sigmoid. This resulted in 72 data forms that were run through the support
vector regressor.

Table 3 shows the best performers in each metric. The best $R^{2}$, 0.051143,
came from the Rochester, Buffalo, Syracuse, and Albany dataset, normalized
with z-score, using the rbf kernel, and with a random state of None. However,
SVR did not perform well on ROC dataset where the invalid values were
calculated in $R^{2}$. Table 3 shows the top 2 best performers in SVR
including the two invalid values in $R^{2}$.

Table 3: Top 2 best metric performers in SVR. Dataset | Normalization | Kernel | Random | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---|---
Mixed | ZScore | rbf | none | 0.051143241 | 0.013033322 | 0.114163574 | 0.314971336
Mixed | MinMax | rbf | 42 | 0.036319941 | 0.011336791 | 0.106474368 | 0.380176974
ROC | None | sigmoid | 42 | invalid | 0.012855096 | 0.113380316 | 0.168736666
ROC | None | sigmoid | 0 | invalid | 0.011151403 | 0.105600204 | 0.171740467

####  4.2.4 DNN Regressor

In total, there were 126 DNN Regressors were trained against the cross of the
6 premade datasets, 7 different amounts of hidden layers, and the 3 Random
States used from before. Seven different hidden layers were: 2, 3, 4, 5, 10,
20, and 30 layers. Each of the hidden layers had the number of inputs equal to
the number of features which for the ROC datasets was 9 and for the
ROC+BUF+SYR+ALB there were 36 inputs at each hidden layer.

Table 4 shows the best two performers in each dataset. The best $R^{2}$,
0.038945, came from mixed dataset, normalized with MinMax, using 5 hidden
layers, and with a random state of None. The best Pearson Correlation metric,
0.267228, came from the mixed dataset, normalized with MinMax, using 3 hidden
layers, and with a random state of None. Even though we used 30 layers deep
neural network, the best performer of $R^{2}$ in ROC dataset was obtained from
the normalization with MinMax, using 4 hidden layers, and a random state of
42.

Table 4: Top 2 best metric performers in DNNr. Dataset | Normalization | Layers | Random | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---|---
Mixed | MinMax | 5 | none | 0.038944969 | 0.010208255 | 0.101035909 | 0.219027928
Mixed | MinMax | 10 | none | 0.0020580 | 0.00821409 | 0.09063162 | 0.05852291
ROC | MinMax | 4 | 42 | 0.036257134 | 0.01133753 | 0.106477838 | 0.24003999
ROC | MinMax | 5 | none | 0.035922961 | 0.010334113 | 0.101656837 | 0.191548407

####  4.2.5 WNN Regressor

Similarly to the WNN Classifier, there were 18 unique WNN Classifiers by
training against the cross of the 6 premade datasets and the 3 Random States
used from before. For each of the measured metrics there were different models
that performed best. Table 5 shows the best performers in each metric. The
best $R^{2}$, 0.132945, came from the Rochester, Buffalo, Syracuse, and Albany
dataset, normalized with ZScore, and with a random state of None. The best
Pearson Correlation metric also came from the mixed dataset, normalized with
MinMax, and with a random state of 42 (not listed on the table). The best
MSE/RMSE values were obtained from the mixed dataset normalized with MinMax,
and a random state of None (not listed on the table). Those values were MSE =
0.007244 and RMSE = 0.085110 with $R^{2}$=0.119940624 , and Pcc=0.357592781.
In ROC dataset, the best $R^{2}$ is obtained with a decent MSE and RMSE.

Table 5: Top 2 best metric performers in WNN regressor. Dataset | Normalization | Random | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---
Mixed | ZScore | none | 0.132945285 | 0.008984219 | 0.094785122 | 0.368640036
Mixed | MinMax | 0 | 0.121693571 | 0.008829052 | 0.093963034 | 0.354665963
ROC | ZScore | 42 | 0.123917918 | 0.010306283 | 0.101519866 | 0.352715523
ROC | MinMax | None | 0.121567081 | 0.009098939 | 0.096595284 | 0.353484533

####  4.2.6 DWNN Regressor

Similarly to the DNN Regressor, 7 different hidden layers were: 2, 3, 4, 5,
10, 20, and 30 layers. Table 6 shows the best performers in each metric. The
best $R^{2}$ for the mixed datasets, 0.181974, came from the mixed dataset,
normalized with ZScore, using 20 hidden layers, and with a random state of 42.
The best $R^{2}$ for the ROC datasets, 0.143873, were obtained from the
Rochester dataset normalized with ZScore, using 20 hidden layers, and a random
state of None.

Table 6: Top 2 best metric performers in DNNr. Dataset | Normalization | Layers | Random | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---|---
Mixed | ZScore | 20 | 42 | 0.181974 | 0.009623 | 0.098098 | 0.426754
Mixed | ZScore | 30 | 42 | 0.181682 | 0.009627 | 0.098116 | 0.428402
ROC | MinMax | 20 | none | 0.143873 | 0.008235 | 0.090748 | 0.379773
ROC | MinMax | 30 | none | 0.108191 | 0.004529 | 0.067294 | 0.361150

####  4.2.7 LSTM Regressor

There were totally 30 LSTM Regressors by training against the cross of the 6
premade datasets, and 5 different sequence lengths. Random states were not
used for LSTM as it depends on chronologically consecutive rows of data. The 5
sequence lengths considered were 3, 5, 7, 9, and 12 previous rows of data.
Table 7 shows the results of the two best performing LSTM regressors over two
datasets. The best $R^{2}$, MSE, and RMSE came from the same model, using the
Rochester dataset, not normalized, and using a sequence length of 5. Those
values were $R^{2}$ = 0.099439, MSE = 0.006424, and RMSE = 0.080149. The best
Pearson Correlation, 0.324997, came from the Rochester dataset, normalized
using ZScore, and using a sequence length of 3. From the mixed dataset, the
regressor with highest $R^{2}$ obtains the lowest MSE and RMSE and the highest
Pcc.

Table 7: Top 2 best metric performers in LSTM. Dataset | Normalization | Sequence | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---
Mixed | MinMax | 12 | 0.093923238 | 0.006489071 | 0.080554771 | 0.310782912
Mixed | None | 12 | 0.072728045 | 0.006640865 | 0.081491504 | 0.272888933
ROC | None | 5 | 0.099438992 | 0.006423882 | 0.080149126 | 0.318826
ROC | ZScore | 12 | 0.088653463 | 0.006526812 | 0.080788686 | 0.313515598

####  4.2.8 LSTM Bi-direction

The Bidirectional Regressor models used the same settings as the LSTM
Regressor models, so 30 unique Bidirectional Regressors were examined. Table 8
shows the results of the two best performing LSTM regressors. The best
$R^{2}$, MSE, RMSE, and Pcc in the Rochester dataset came from the same
regressor, normalized using MinMax, and using a sequence length of 12. Those
values were $R^{2}$ = 0.092928, MSE = 0.006496, and RMSE = 0.080599. The best
$R^{2}$, MSE, RMSE, and Pcc in the mixed dataset also came from the same
regressor that was normalized using MinMax, and using a sequence length of 3.

Table 8: Top 2 best metric performers in LSTM Bi-direction. Dataset | Normalization | Sequence | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---
Mixed | MinMax | 3 | 0.067291762 | 0.006645372 | 0.081519154 | 0.29581004
Mixed | MinMax | 7 | 0.06109682 | 0.006705727 | 0.081888505 | 0.28511238
ROC | MinMax | 12 | 0.09292772 | 0.006496201 | 0.080599012 | 0.320602646
ROC | None | 12 | 0.085119002 | 0.006552125 | 0.080945195 | 0.311874173

####  4.2.9 GRU Regressor

The GRU Regressor models used the same settings as the LSTM and Bidirectional
Regressor models, so 30 unique GRU Regressors. The best metrics in ROC dataset
were normalized using ZScore, and a sequence length of 3. Those values were
$R2$ = 0.088520, Pearson Correlation = 0.323445, MSE = 0.006494, and RMSE =
0.080586. The best metrics in mixed dataset were normalized using MinMax, and
a sequence length of 9. Those values were $R2$ = 0.073344684, Pearson
Correlation = 0.293189287, MSE = 0.006625979, and RMSE = 0.081400116. Table 9
shows the top two performing GRU Regressor Models and their respective
performance metrics.

Table 9: Top 2 best metric performers in GRU. Dataset | Normalization | Sequence | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---
Mixed | MinMax | 9 | 0.073344684 | 0.006625979 | 0.081400116 | 0.293189287
Mixed | MinMax | 5 | 0.067207667 | 0.006653795 | 0.081570797 | 0.29401964
ROC | ZScore | 3 | 0.088519946 | 0.006494126 | 0.080586137 | 0.323445285
ROC | MinMax | 5 | 0.079144128 | 0.00656865 | 0.081047207 | 0.301406376

####  4.2.10  Overall Regression Results

Table 10 demonstrates the best performer of each models in the mixed
ROC+BUF+SYR+ALB dataset. It shows that KNN gives the best $R^{2}$ that means
that the prediction model best fits the observations at the same time its MSE
and RMSE are superior, similar to those metrics in the models of LSTM, GRU,
and Bidirectional LSTM, 0.006806884 vs (0.006489071, 0.006625979, and
0.006645372). Also, KNN’s Pcc shows the highest value among the models, which
means a high correlation coefficient between the prediction model and the
dataset. The $R^{2}$ value of KNN outperforms the rest models. The second
place is DWNNr with 20 layers and ZScore normalization.

Recurrent neural network models were expected to have a better performance in
sequential data prediction. However, recurrent neural network models such as
LSTM, GRU, and bidirectional LSTM have not shown a strong indicator in $R^{2}$
and Pcc. The experimental results partially match the previous researcher’s
findings. [Liu et al. (2020)] revealed that the rainfall time series have very
short term memory characteristics. It sheds some lights on why LSTM models
does not give the best performance in the precipitation forecasting.

Table 10: Overall Performance in ROC+BUF+SYR+ALB dataset Model | Normalization | Random | Parameter | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---|---
KNN | ZScore | None | 42 nodes | 0.209561502 | 0.006806884 | 0.082503842 | 0.464431631
DWNNr | ZScore | 42 | 20 layers | 0.181974116 | 0.009623306 | 0.098098453 | 0.426754455
Linear | ZScore | None | N/A | 0.1464424 | 0.007025658 | 0.083819199 | 0.388019644
WNNr | ZScore | None | N/A | 0.132945285 | 0.008984219 | 0.094785122 | 0.368640036
LSTM | MinMax | N/A | 12 units | 0.093923238 | 0.006489071 | 0.080554771 | 0.310782912
GRU | MinMax | N/A | 9 units | 0.073344684 | 0.006625979 | 0.081400116 | 0.293189287
Bidirect | MinMax | N/A | 3 units | 0.067291762 | 0.006645372 | 0.081519154 | 0.29581004
DNNr | MinMax | None | 5 layers | 0.038944969 | 0.010208255 | 0.101035909 | 0.219027928
SVR | ZScore | 42 | rbf kernel | 0.036319941 | 0.011336791 | 0.106474368 | 0.380176974

In the ROC dataset, KNN is ranked the second place with a little lower value
in $R^{2}$ than that in the linear regression model, showing the performance
0.178324011 vs. 0.168353653. However, KNN is much better in MSE and RMSE than
those of the first place. Comparing with the third place DWNNr, all metrics of
KNN are superior. The ROC dataset is not so big as the mixed dataset that KNN
may be a little inferior in $R^{2}$ and Pcc than the first place but in other
metrics KNN exhibited its strong prediction ability in rainfall forecasting.

Table 11: Overall Performance in ROC dataset Model | Normalization | Random | Parameter | $R^{2}$ | MSE | RMSE | Pcc
---|---|---|---|---|---|---|---
Linear | None | None | N/A | 0.178324011 | 0.006913207 | 0.083145696 | 0.422665921
KNN | MinMax | None | 51 nodes | 0.168353653 | 0.005346776 | 0.07312165 | 0.411062473
DWNNr | ZScore | None | 20 layers | 0.143873162 | 0.008235172 | 0.090747846 | 0.379772590
WNNr | ZScore | 42 | N/A | 0.123917918 | 0.010306283 | 0.101519866 | 0.352715523
LSTM | None | N/A | 5 units | 0.099438992 | 0.006423882 | 0.080149126 | 0.318826
Bidirect | MinMax | N/A | 12 units | 0.09292772 | 0.006496201 | 0.080599012 | 0.320602646
GRU | ZScore | N/A | 3 units | 0.088519946 | 0.006494126 | 0.080586137 | 0.323445285
DNNr | MinMax | 42 | 4 layers | 0.036257134 | 0.01133753 | 0.106477838 | 0.24003999
SVR | None | 42 | sigmoid | invalid | 0.012855096 | 0.113380316 | 0.168736666

##  5 Conclusion and Discussion

Rainfall forecasting plays an important role in our daily lives, especially
agriculture and related activities. We have designed an integrated tool of
machine learning modeling for regional rainfall forecasting. It applied the
state-of-the-art machine learning algorithms including Deep Neural Network
(DNN), Reservoir Computing (RC), Wide Neural Network (WNN), Deep and Wide
Neural Network (DWNN), Long Short Term Memory (LSTM), Support Vector Machine
(SVM), and K-Nearest Neighbor (KNN). Not only the classification models but
also the regression models were implemented. In the classification models, we
used the classification models to forecast if the precipitation will occur in
next hour. In the regression models, we used the regression models to predict
the regression value of the rainfall. Also, we adopted two normalization
methods MinMax and ZScore and compared their performance with those of non-
normalized models. The results shows that KNN model with ZScore normalization
can outperform over other models according to the metrics of $R^{2}$, MSE,
RMSE, and Pcc.

Time series prediction models such as LSTM, GRU, and LSTM-bidirection have not
shown satisfactory performance in rainfall forecasting. On the contrary, the
linear based models such as DWNN, linear models, and WNN have much better
performance than complicated models such as DNN, SVR, and LSTM. The reason
behind the ranking is that the weather data unlike other data contains many
uncertain outliers, errors, and missings that may fuel lots of distractions to
the precise models such as neural networks, and time series models. Therefore,
rainfall forecasting requires more robust models to analyze those weather
datasets. Fortunately, through the experimental results, we can find that the
KNN model, a model that has been underestimated, is robust enough to be
qualified for this type of job. In the future, we will focus on refining the
rainfall forecasting models, especially those methods that can detect the
anomaly in order to resist their impact on weather datasets.

## References

  * Bajpai and Bansal (2021) Bajpai, V., and A. Bansal. 2021.  “A Deep and Wide Neural Network-based Model for Rajasthan Summer Monsoon Rainfall (RSMR) Prediction”.  arXiv preprint arXiv:2103.02157.
  * Bojang et al. (2020) Bojang, P. O., T.-C. Yang, Q. B. Pham, and P.-S. Yu. 2020.  “Linking singular spectrum analysis and machine learning for monthly rainfall forecasting”.  Applied Sciences vol. 10 (9), pp. 3224.
  * Canchala et al. (2020) Canchala, T., W. Alfonso-Morales, Y. Carvajal-Escobar, W. L. Cerón, and E. Caicedo-Bravo. 2020.  “Monthly rainfall anomalies forecasting for southwestern Colombia using artificial neural networks approaches”.  Water vol. 12 (9), pp. 2628.
  * Cheng et al. (2016) Cheng, H.-T., L. Koc, J. Harmsen, T. Shaked, T. Chandra, H. Aradhye, G. Anderson, G. Corrado, W. Chai, M. Ispir et al. 2016.  “Wide & deep learning for recommender systems”.  In Proceedings of the 1st workshop on deep learning for recommender systems, pp. 7–10.
  * Chhetri et al. (2020) Chhetri, M., S. Kumar, P. Pratim Roy, and B.-G. Kim. 2020.  “Deep BLSTM-GRU Model for Monthly Rainfall Prediction: A Case Study of Simtokha, Bhutan”.  Remote Sensing vol. 12 (19), pp. 3174.
  * Cho et al. (2014) Cho, K., B. Van Merriënboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio. 2014.  “Learning phrase representations using RNN encoder-decoder for statistical machine translation”.  arXiv preprint arXiv:1406.1078.
  * Coble and Yu (2020) Coble, N. J., and N. Yu. 2020.  “A Reservoir Computing Scheme for Multi-class Classification”.  In Proceedings of the 2020 ACM Southeast Conference, pp. 87–93.
  * French et al. (1992) French, M. N., W. F. Krajewski, and R. R. Cuykendall. 1992.  “Rainfall forecasting in space and time using a neural network”.  Journal of Hydrology vol. 137 (1), pp. 1–31.
  * Gallicchio and Micheli (2017) Gallicchio, C., and A. Micheli. 2017.  “Deep echo state network (deepesn): A brief survey”.  arXiv preprint arXiv:1712.04323.
  * Hasan et al. (2015) Hasan, N., N. C. Nath, and R. I. Rasel. 2015.  “A support vector regression model for forecasting rainfall”.  In 2015 2nd International Conference on Electrical Information and Communication Technologies (EICT), pp. 554–559.  IEEE.
  * Hochreiter and Schmidhuber (1997) Hochreiter, S., and J. Schmidhuber. 1997.  “Long short-term memory”.  Neural computation vol. 9 (8), pp. 1735–1780.
  * Hu et al. (2013) Hu, J., J. Liu, Y. Liu, and C. Gao. 2013.  “EMD-KNN model for annual average rainfall forecasting”.  Journal of Hydrologic Engineering vol. 18 (11), pp. 1450–1457.
  * Huang et al. (2017) Huang, M., R. Lin, S. Huang, and T. Xing. 2017.  “A novel approach for precipitation forecast via improved K-nearest neighbor algorithm”.  Advanced Engineering Informatics vol. 33, pp. 89–95.
  * Hung et al. (2009) Hung, N. Q., M. S. Babel, S. Weesakul, and N. Tripathi. 2009.  “An artificial neural network model for rainfall forecasting in Bangkok, Thailand”.  Hydrology and Earth System Sciences vol. 13 (8), pp. 1413–1425.
  * Jan et al. (2008) Jan, Z., M. Abrar, S. Bashir, and A. M. Mirza. 2008.  “Seasonal to inter-annual climate prediction using data mining KNN technique”.  In International Multi Topic Conference, pp. 40–51.  Springer.
  * Koizumi (1999) Koizumi, K. 01 Feb. 1999.  “An Objective Method to Modify Numerical Model Forecasts with Newly Given Weather Data Using an Artificial Neural Network”.  Weather and Forecasting vol. 14 (1), pp. 109 – 118.
  * Liu et al. (2020) Liu, M., Y. Huang, Z. Li, B. Tong, Z. Liu, M. Sun, F. Jiang, and H. Zhang. 2020\.  “The applicability of LSTM-KNN model for real-time flood forecasting in different climate zones in China”.  Water vol. 12 (2), pp. 440.
  * Luk et al. (2001) Luk, K. C., J. Ball, and A. Sharma. 2001.  “An application of artificial neural networks for rainfall forecasting”.  Mathematical and Computer Modelling vol. 33 (6), pp. 683–693.
  * Pölsterl et al. (2019) Pölsterl, S., I. Sarasua, B. Gutiérrez-Becker, and C. Wachinger. 2019.  “A wide and deep neural network for survival analysis from anatomical shape and tabular clinical data”.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pp. 453–464.  Springer.
  * Todaro (2018) Todaro, M. 2018.  “Assessing Irrigation in Western New York State”.  Thesis SUNY ESF.
  * Toth et al. (2000) Toth, E., A. Brath, and A. Montanari. 2000.  “Comparison of short-term rainfall prediction models for real-time flood forecasting”.  Journal of hydrology vol. 239 (1-4), pp. 132–147.
  * Vlachas et al. (2020) Vlachas, P. R., J. Pathak, B. R. Hunt, T. P. Sapsis, M. Girvan, E. Ott, and P. Koumoutsakos. 2020.  “Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics”.  Neural Networks vol. 126, pp. 191–217.
  * Wu (2009) Wu, J. 2009.  “A Novel Artificial Neural Network Ensemble Model Based on K–Nearest Neighbor Nonparametric Estimation of Regression Function and Its Application for Rainfall Forecasting”.  In 2009 international joint conference on computational sciences and optimization, Volume 2, pp. 44–48.  IEEE.
  * Yang et al. (2020) Yang, M., H. Wang, Y. Jiang, X. Lu, Z. Xu, and G. Sun. 2020.  “Geca proposed ensemble–knn method for improved monthly runoff forecasting”.  Water Resources Management vol. 34 (2), pp. 849–863.
  * Yen et al. (2019) Yen, M.-H., D.-W. Liu, Y.-C. Hsin, C.-E. Lin, and C.-C. Chen. 2019.  “Application of the deep learning for the prediction of rainfall in Southern Taiwan”.  Scientific reports vol. 9 (1), pp. 1–9.


# Improving Energy Efficiency of Massive MIMO Relay Systems using Power
Bisection Allocation for Cell-Edge Users

Mohammad Koolivand†, Mohammad Hossein Bahonar†, Mohammad Sadegh Fazel†
† Deptartment of Electrical and Computer Engineering, Isfahan University of
Technology, Isfahan, Iran 84156-83111
Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper studies downlink energy efficiency (EE) of a cellular massive
Multiple Input Multiple Output (MIMO) system with massive number of antennas
at both base station (BS) and relay considering cell-edge users which is less
investigated in the literature. Some improving methods for base station (BS)
and relay power allocation under quality of service constraints are proposed.
The EE optimization problem is provided over BS and relay power variables
(power dimensions). Considering the quasi-concavity property of the EE
function relative to each power variable, we propose to use a power bisection
algorithm (PBA) in one dimension, followed by the exhaustive search in the
other dimension (ODS algorithm) which is called PB-ODS algorithm and also a
method to limit the range of power variables. The results show that the
performance of the PB-ODS algorithm approaches the performance of optimal
solution which is exhaustive search in both dimensions while it has a lower
complexity. In addition to that we suggest to use the PBA in two dimensions
alternatively as the sub-optimal alternative optimization (AOP) algorithm. The
complexity can highly be reduced using the AOP algorithm with a slight but
acceptable degradation in performance which makes the algorithm much suitable
for practical use.

###### Index Terms:

Energy efficiency, bisection algorithm, massive MIMO relay system, alternative
algorithm, quasi-concave function.

##  I Introduction

Due to the ever increasing rate demand in cellular netwroks, new technologies
and architectures such as heterogeneous networks[1], device-to-device
communication[2], and Massive MIMO[3] have been proposed in order to increase
the overall sum-rate or spectral efficiency (SE) of wireless systems. Massive
MIMO along with cooperative relaying systems, which is known as one of the
most important technologies for increasing network capacity and network
coverage, can greatly help to meet some future network’s requirements and can
be considered as a good choice for the next generation of wireless
communication systems. In the early researches[4, 5] , some challenges such as
system capacity have been addressed.

In [6], the sum rate of a system where several pairs of single-antenna users
communicate with each other through a one-way amplify and forward (AF) MIMO
relay is obtained. In [7, 8], considering a two-way AF MIMO relay system,
maximizing the total network rate is investigated. In [9] a full-duplex model
for a one-way AF relay system with maximum ratio combining/transmission
(MRC/MRT) is proposed and a lower bound for the total transmission rate of the
network is obtained using Jensen inequality. The scenario considered in [3] is
a one-way massive MIMO AF relay system in which its energy efficiency (EE) is
investigated and some methods based on Random Matrix Theory are proposed to
improve it. In [10], the achievable spectral efficiency (SE) of a multi-cell
AF multi-way massive MIMO relay network in presence of some practical
transmission impairments, including imperfect channel estimation or outdated
channel state information (CSI) has been investigated.

In recent years, some researchers have been conducting investigations on
massive MIMO relay systems in cellular networks. Researchers in [11] have
considered a cellular network with MIMO relay and investigated its total SE.
Different configuretions for single-cell massive MIMO relay systems have been
performed recently in [12]. In [13], authors investigate the tradeoff between
SE and EE in AF relay networks. The system model of [14] includes a large
number of single-antenna users, which communicate with a massive MIMO base
station (BS) with aid of single-antenna AF relays. The EE and SE of the system
is obtained using the bisection algorithm and is compared with the exhaustive
search method.

Due to increasing power consumption in recent years, design of energy
efficient systems have attracted a lot of researches. Because of the
importance of EE in massive MIMO systems, we have investigated the combination
of massive MIMO and cooperative relaying systems with energy efficient
structures for 5G and beyond cellular networks. Few papers have investigated
the EE in a cellular massive MIMO relay-aided network while the SE and the
capacity of the system is usually examined; for example [15] investigates
massive MIMO two-way relays in cellular networks. So in this paper we study
the EE of a single-cell massive MIMO relay system and model the power
allocation problem using a non-convex optimization problem. We aim to serve
cell-edge users in a downlink transmission scenario. Circuit power consumption
is considered in the power consumption model. Due to the non-convex form of
the optimization problem, its optimal solution is 2-D exhaustive search (TDS)
algorithm over BS and relay power variables but it is a highly computationally
complex approach. It is shown that our EE function is a quasi-concave function
with respect to each power variable, therefore we propose a method to limit
the range of power variables and also PB-ODS method which uses a power
bisection algorithm (PBA) with respect to one of the power variables following
an exhaustive search with respect to the other power variable. It is also
suggested to use the PBA with respect to two power variables alternatively as
the sub-optimal alternative optimization (AOP) algorithm. This method has a
slight but acceptable degradation in performance while its computational
complexity has dramatically decreased.

The rest of the paper is organized as follows. The system model for a massive
MIMO relay cellular system and some signaling parameters are introduced in
section II. Section III describes the EE optimization problem formulation.
Using some features of the objective function, we propose an iterative
algorithm for solving the problem in section IV. The simulation results are
presented in section V and finally, section VI concludes the paper.

Notations\- Throughout , matrices are written in capital bold letters such as
$\mathbf{A}$, vectors are in bold lower case such as $\mathbf{a}$, and scalars
as plain letter such as $A$ or $a$. $\mathbf{A}^{T}$ denotes the transpose of
$\mathbf{A}$, $\mathbf{A}^{h}$ is the conjugate transpose of $\mathbf{A}$,
$\mathbf{A}^{-1}$ represents the inverse of $\mathbf{A}$,
$\mathbf{A}^{\dagger}$ corresponds to the pseudo- inverse of $\mathbf{A}$, and
$tr(\mathbf{A})$ denotes its trace. $\mathbb{E}(.)$ represents the expectation
operator.

##  II System model

Consider a single cell containing a BS at its center with massive number of
antennas $N_{b}$, a MIMO AF relay with $N_{r}$ antennas, and $K$ single-
antenna users that are randomly distributed in a cell according to Fig. 1. The
number of relay antennas is lower than the number of BS antennas
$(N_{b}=\alpha N_{r}\leavevmode\nobreak\ ,\leavevmode\nobreak\ \alpha\geq 1)$.
The relay operates in half-duplex mode which means that data transmission from
the BS to its corresponding users is performed in two consecutive time slots.
Achievable data rate and power consumption model are presented in the
following.

Figure 1: Downlink transmission through a massive MIMO relay in a single cell
system.

###  II-A Achiveable Sum Rate

According to Fig. 1, in the first time slot the BS sends data to the relay.
The data at the BS is denoted as
$\mathbf{\tilde{s}}_{b}\in\mathcal{C}^{K\times 1}$ such that
$\mathbb{E}(\mathbf{\tilde{s}}_{b}\mathbf{\tilde{s}}_{b}^{h})=\mathbf{I}_{K}$,
so the signal that is sent from the BS can be expressed as

$\mathbf{x}_{b}=\mathbf{F}_{b}\mathbf{s}_{b},$ (1)

where $\mathbf{F}_{b}\in\mathcal{C}^{N_{b}\times N_{r}}$ is the precoding
matrix at the BS. The matrix $\mathbf{F}_{b}$ for ZF proccessing is equal to
$\alpha_{b}\mathbf{H}_{BR}^{\dagger}$ where
$\mathbf{H}_{BR}=\sqrt{d_{BR}}{\mathbf{G}}_{BR}$ is the channel matrix between
the BS and relay. $d_{BR}$ and $\mathbf{G}_{BR}$ are the large-scale log-
normal fading coefficient and the small-scale zero mean unit variance complex
normal fading matrix, respectively. $\mathbf{s}_{b}$ is formed by adding a
$(N_{r}-K)\times 1$ dimension zero vector to $\tilde{\mathbf{s}}_{b}$, so we
have $\mathbf{s}_{b}=\left(\begin{smallmatrix}\tilde{\mathbf{s}}_{b}\\\
0_{(N_{r}-K)\times 1}\end{smallmatrix}\right)$ such that
$\mathbb{E}{(\mathbf{s}_{b}{\mathbf{s}_{b}}^{h})}=\left(\begin{smallmatrix}\mathbf{I}_{K}&0\\\
0&0\end{smallmatrix}\right)$. It is also assumed that full channel state
information (CSI) is known at the BS. The power normalization factor
($\alpha_{b}$) is derived from the fact that the instantaneous transmission
power of the BS per timeslot is equal to $P_{b}$ and can be expressed as

$\displaystyle\begin{split}tr[\mathbb{E}(\mathbf{x}_{b}\mathbf{x}_{b}^{h})]=P_{b}&\Rightarrow
tr[(\alpha_{b}\mathbf{H}_{BR}^{\dagger})\mathbb{E}(\mathbf{s}_{b}\mathbf{s}_{b}^{h})(\alpha_{b}\mathbf{H}_{BR}^{\dagger})^{h}]=P_{b}\\\
&\alpha_{b}^{2}tr[\tilde{\mathbf{H}}_{BR}^{\dagger}(\tilde{\mathbf{H}}_{BR}^{\dagger})^{h}]=P_{b}\\\
&\alpha_{b}=\sqrt{\frac{P_{b}}{tr(\tilde{\mathbf{H}}_{BR}^{\dagger}(\tilde{\mathbf{H}}_{BR}^{\dagger})^{h})}}.\end{split}$
(2)

Since the $\mathbf{s}_{b}$ vector contains $(N_{r}-K)$ zero elements, the
matrix $\tilde{\mathbf{H}}_{BR}^{\dagger}$ is derived by selecting the first
$K$ columns of the matrix $\mathbf{H}_{BR}^{\dagger}$. The normalization
factor is calculated by considering the vector $\tilde{\mathbf{s}}_{b}$ and
the matrix $\tilde{\mathbf{H}}_{BR}^{\dagger}$. The transmitted data from the
BS is received at the relay after passing through the BS-to-relay fading
channel. So the signal vector received at the relay can be expressed as

$\displaystyle\mathbf{y}_{r}=\mathbf{H}_{BR}\mathbf{F}_{b}\mathbf{s}_{b}+\mathbf{z}_{r}=\alpha_{b}\mathbf{s}_{b}+\mathbf{z}_{r}.$
(3)

It is assumed that the distance between the BS antennas and relay antennas is
much greater than the distance between the relay antennas with each other (and
also the distance between the BS antennas from each other). Also,
$\mathbf{z}_{r}\sim\mathcal{CN}(0,\sigma^{2}_{r}\mathbf{I}_{N_{r}})$ is the
additive white Gaussian (AWG) noise. The relay multiplies $\mathbf{y}_{r}$ by
the beamforming matrix $\mathbf{W}_{r}\in\mathcal{C}^{N_{r}\times K}$ for
amplification and ZF processing. Therefore, the data sent from the relay is
equal to $\mathbf{x}_{r}=\mathbf{W}_{r}\tilde{\mathbf{y}}_{r}$ where the
vector $\tilde{\mathbf{y}}_{r}$ is obtained by choosing the first $K$ elements
of vector $\mathbf{y}_{r}$. The processing matrix is equal to
$\mathbf{W}_{r}=\alpha_{r}\mathbf{H}_{RU}^{\dagger}$, where $\mathbf{H}_{RU}$
is the channel matrix between the relay and users and
$\mathbf{H}_{RU}=\mathbf{G}_{RU}\mathbf{D}^{1/2}$. $\mathbf{D}$ is a diagonal
matrix containing large-scale log-normal fading coefficients and
$\mathbf{G}_{RU}$ is a small-scale zero mean unit variance complex normal
fading channel matrix. The signal sent from the relay is equal to

$\begin{split}\mathbf{x}_{r}&=\mathbf{W}_{r}\tilde{\mathbf{y}}_{r}\\\
&=\alpha_{b}\alpha_{r}\mathbf{H}_{RU}^{\dagger}\tilde{\mathbf{s}}_{b}+\alpha_{r}\mathbf{H}_{RU}^{\dagger}\tilde{\mathbf{z}}_{r},\end{split}$
(4)

where the vector $\tilde{\mathbf{z}}_{r}$ is obtained by choosing the first
$K$ elements of vector $\mathbf{z}_{r}$. The normalization factor $\alpha_{r}$
is chosen such that $tr[\mathbb{E}(\mathbf{x}_{r}\mathbf{x}_{r}^{h})]=P_{r}$
and can be expressed as

$\alpha_{r}=\sqrt{\frac{P_{r}}{(\alpha_{b}^{2}+\sigma^{2}_{r})tr(\mathbf{H}_{RU}\mathbf{H}_{RU}^{h})^{-1}}}.$
(5)

After all, the received signal for the $k$th user is

${y}_{u,k}=\alpha_{b}\alpha_{r}\tilde{{s}}_{b,k}+\alpha_{r}\tilde{{z}}_{r,k}+{z}_{u,k},$
(6)

where ${z}_{u,k}$ is zero mean, $\sigma^{2}_{u}$ variance AWG noise at users’
antennas. Therefore, the instantaneous received signal to interference plus
noise ratio (SINR) of the $k$th user is

$\gamma_{u,k}=\frac{\alpha_{b}^{2}\alpha_{r}^{2}}{\alpha_{r}^{2}\sigma_{r}^{2}+\sigma_{u}^{2}}.$
(7)

Note that the equation (7) holds if users know the values of $\alpha_{b}$ and
$\alpha_{r}$. Finally, the overall sum-rate for the downlink transmission can
be calculated as follows

$\displaystyle R_{sum}$ $\displaystyle=\sum_{k=1}^{K}R_{k}$ (8) $\displaystyle
R_{k}$ $\displaystyle=\frac{W}{2}log_{2}(1+\gamma_{u,k}),$ (9)

where the factor $1/2$ is due to the half-duplex relay transmission mode and
$W$ is the transmission bandwidth.

###  II-B Power Consumption

Total power consumption model can be described as,

$P_{total}=\zeta_{b}P_{b}+\zeta_{r}P_{r}+P_{sys},$ (10)

where $\zeta_{b}$ and $\zeta_{r}$ are the BS and relay amplifier efficiencies,
respectively. $P_{sys}$ is the system power including power consumption of the
RF circuits, digital processors, coding, etc. This power can be divided to the
following parts:

  * •

RF circuit power consumption: The power consumption of the downlink
transmission consists of two parts: BS-to-relay and relay-to-user circuit
power consumption denoted by $P_{c1}=N_{b}P_{tx}+N_{r}P_{rx}$ and
$P_{c2}=N_{r}P_{tx}+KP_{rx}$, respectively. $P_{tx}$ and $P_{rx}$ are RF
circuit power consumption for each antenna in the transmitter and receiver,
respectively. Due to the half-duplex relay transmission mode, the total RF
power consumption can be expressed as

$P_{c}=\frac{1}{2}(P_{c1})+\frac{1}{2}(P_{c2}).$ (11)

  * •

Coding/decoding power: Similarly, the power consumption for coding/decoding
procedure can be expressed as

$P_{cd}=\frac{K}{2}(P_{cod}+P_{dec}),$ (12)

where $P_{cod}$ and $P_{dec}$ correspond to coding and decoding power
consumptions for each user, respectively.

  * •

Precoding power consumption: Similar to [16], the power consumption associated
with the precoding procedure at the BS and relay are denoted by $P_{pc,BS}$
and $P_{pc,R}$, respectively and can be expressed as

$P_{pc,BS}=\frac{3N_{b}^{2}N_{r}+2N_{b}N_{r}}{2LT}+\frac{N_{r}^{3}}{3LT},$
(13) $P_{pc,R}=\frac{3N_{r}^{2}K+2N_{r}K}{2LT}+\frac{N_{r}^{3}}{3LT},$ (14)

where $T$ is the coherence time of the channel. $L$ is BS operation capability
per joule. Therefore the total precoding power consumption is
$P_{pc}=1/2(P_{pc,BS}+P_{pc,R})$. So the total system power consumption can be
calculated from the following equation

$P_{sys}=P_{c}+P_{cd}+P_{pc}+P_{0},$ (15)

where $P_{0}$ is a constant power consumption related to backhaul traffic and
cooling.

##  III Problem Formulation

According to the sum-rate equation (8) and the power consumption model (15) ,
the EE of the downlink transmission can be expressed as

$\eta_{EE}=\frac{R_{sum}}{P_{total}}=\frac{\frac{W}{2}\sum_{k=1}^{K}log_{2}(1+\gamma_{u,k})}{\zeta_{b}P_{b}+\zeta_{r}P_{r}+P_{sys}}.$
(16)

The EE value represents the number of transmitted bits for the consumption of
one joule of energy and our goal is to find optimal power allocation for BS
and relay in order to maximize the EE function (16) considering the minimum
and maximum transmission power constraints of the BS and relay alongside with
the QoS constraints of users. For each user, its QoS constraint is represented
by its minimum achievable data rate, $R_{min,k}$. Therefore the EE
optimization problem can be formulated as

$\displaystyle\max_{P_{b},P_{r}}\leavevmode\nobreak\ $
$\displaystyle\frac{\frac{W}{2}\sum_{k=1}^{K}log_{2}(1+\gamma_{u,k})}{\zeta_{b}P_{b}+\zeta_{r}P_{r}+P_{sys}}$
(17a) $\displaystyle s.t.$ $\displaystyle 0\leq P_{b}\leq P_{b_{max}},$ (17b)
$\displaystyle 0\leq P_{r}\leq P_{r_{max}},$ (17c) $\displaystyle R_{k}\geq
R_{min,k},\leavevmode\nobreak\ k=1,...,K$ (17d)

where $P_{b_{max}}$ and $P_{r_{max}}$ correspond to the maximum transmission
powers of the BS and relay, respectively. Constraints (17b) and (17c) indicate
the minimum and maximum transmission powers of the BS ans relay, respectively.
Constraint (17d) indicates the QoS constraints of users.

##  IV Proposed Method

By solving the EE optimization problem in (17a), the optimal values of BS and
relay transmission powers can be found. The ideal but extremely complex
approach is the exhaustive search method. This optimal method searches among
all possible values of transmission powers but it may not be a good practical
solution due to its very high complexity. Another solution is to optimize the
relay and BS power according to behavior of the EE function. The structure of
our EE function is similar to the EE function of [13]. It is stated in [13]
that their EE function is a quasi-concave function. Therefore, it is
reasonable to conclude that our EE function is also a quasi-concave function
relative to the relay or BS power variables. Simulation results also confirm
this conclusion. Using this property, the optimum value of the EE function can
be found using a power bisection algorithm. So we propose the PBA algorithm
for relay or BS power allocation based on the quasi-concavity property of the
EE function in order to maximize the EE function.

1:$\rm{P_{r_{max}}\leavevmode\nobreak\
,P_{b}^{(0)},\epsilon_{r}>0,\kappa,P_{r_{min}}}$

2:$\rm{\eta_{EE}^{*},\leavevmode\nobreak\ P_{r}^{*}}$

3:$\rm{initialize\leavevmode\nobreak\ P^{(1)}=P_{r_{min}}}$

4:$\rm{calculate\leavevmode\nobreak\ \leavevmode\nobreak\
\eta_{EE}\leavevmode\nobreak\ (for\leavevmode\nobreak\
P_{r}=P^{(1)}\leavevmode\nobreak\ \&\leavevmode\nobreak\
P_{b}^{(0)})\leavevmode\nobreak\ from\leavevmode\nobreak\
eq\leavevmode\nobreak\ \eqref{eq_ee}}$

5:if $\rm{0\geq\frac{\partial\eta_{EE}}{\partial P_{r}}}$ then

6: $\rm{P_{r}^{*}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
P^{(1)}}$

7: $\rm{\eta_{EE}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
\eta_{EE}\leavevmode\nobreak\ (for\leavevmode\nobreak\ P_{r}=P^{(1)})}$

8: break

9:else

10: while $\rm{\leavevmode\nobreak\ 0<\leavevmode\nobreak\
\frac{\partial\eta_{EE}}{\partial P_{r}}\leavevmode\nobreak\
\leavevmode\nobreak\ \&\&\leavevmode\nobreak\ \leavevmode\nobreak\
P^{(1)}<P_{r_{max}}}$ do

11: $\rm{P^{(2)}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
P^{(1)}}$

12: $\rm{P^{(1)}\leavevmode\nobreak\ \longleftarrow min(\kappa
P^{(1)},P_{r_{max}}),\leavevmode\nobreak\ \kappa>1}$

13: $\rm{calculate\leavevmode\nobreak\ \eta_{EE}^{1}\leavevmode\nobreak\
for\leavevmode\nobreak\ P_{r}=P^{(1)},\leavevmode\nobreak\
P_{b}^{(0)}\leavevmode\nobreak\ from\leavevmode\nobreak\ \eqref{eq_ee}}$

14: end while

15: if $\rm{0\leq\frac{\partial\eta_{EE}}{\partial P_{r}}}$ then

16: $\rm{P_{r}^{*}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
P^{(1)}}$

17: $\rm{\eta_{EE}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
\eta_{EE}\leavevmode\nobreak\ (for\leavevmode\nobreak\ P_{r}=P^{(1)})}$

18: break

19: end if

20:end if

21:$\rm{calculate\leavevmode\nobreak\ \eta_{EE}^{2}\leavevmode\nobreak\
for\leavevmode\nobreak\ P_{r}=P^{(2)},P_{b}^{(0)}\leavevmode\nobreak\
from\leavevmode\nobreak\ \eqref{eq_ee}}$

22:while $\rm{|\eta_{EE}^{2}-\eta_{EE}^{1}|>\epsilon\leavevmode\nobreak\
\leavevmode\nobreak\ }$ do

23: $\rm{\overline{P}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
\frac{P^{(1)}+P^{(2)}}{2}}$

24: $\rm{calculate\leavevmode\nobreak\ \eta_{EE}\leavevmode\nobreak\
for\leavevmode\nobreak\ P_{r}=\overline{P},P_{b}^{(0)}\leavevmode\nobreak\
from\leavevmode\nobreak\ \eqref{eq_ee}}$

25: if $\rm{0\leq\frac{\partial\eta_{EE}}{\partial P_{r}}}$ then

26: $\rm{P^{(2)}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
\overline{P}}$

27: $\rm{\eta_{EE}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
\eta_{EE}^{2}}$

28: else

29: $\rm{P^{(1)}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
\overline{P}}$

30: $\rm{\eta_{EE}\leavevmode\nobreak\ \longleftarrow\leavevmode\nobreak\
\eta_{EE}^{1}}$

31: end if

32:end while

33:$\rm{Output\leavevmode\nobreak\ \leavevmode\nobreak\ :\leavevmode\nobreak\
\leavevmode\nobreak\ P_{r}^{*}=\overline{P}\leavevmode\nobreak\
\leavevmode\nobreak\ ,\leavevmode\nobreak\ \leavevmode\nobreak\
\eta_{EE}^{*}=\eta_{EE}}$

Algorithm 1 PBA for relay Power

Unfortunately the EE function is not a quasi-concave function of both BS and
relay power simultaneously[13]. So the PBA should be used only for one power
variable (power dimension) at each time. So it is proposed to apply the PBA in
one dimension along with the 1-D exhaustive search in the other dimension
which is called the 1-D Search (ODS) algorithm. Therefore there can be two
types of PB-ODS (PBR-ODSB and PBB-ODSR) algorithms depending on the power
variable which PBA is applied to. PBR-ODSB (PBB-ODSR) correspond to the
algorithm where PBA is applied to relay (BS) power and ODS is applied to BS
(realy) power.

####  IV-1 Relay power optimization

When the BS power is assumed to be constant in equation (17a) the relay power
optimization can be defined as

$\displaystyle\max_{P_{r}}\frac{\frac{W}{2}\sum_{k=1}^{K}log_{2}(1+\gamma_{u,k})}{\zeta_{b}P_{b}+\zeta_{r}P_{r}+P_{sys}}$
$\displaystyle s.t.\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ 0\leq P_{r}\leq P_{r_{max}},$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ R_{k}\geq
R_{min,k},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ k=1,...,K.$ (18)

Since the objective function is quasi-concave with respect to the relay power,
the PBA can be used to find its optimum value as described in algorithm 1
where $\kappa$ is the convergence step and $\rm{P_{r_{min}}}$ is the lowest
feasible relay transmission power which is used for faster convergence of the
algorithm by limiting the relay transmission power variable. Using the QoS
inequality constraints of each user and some mathematical manipulations,
$\rm{P_{r_{min}}}$ can be expressed as

$P_{r_{min}}=\max_{k=1,...,K}\frac{(2^{2R_{min,k}}-1)M\sigma^{2}_{u}}{\alpha_{b}^{2}-(2^{2R_{min,k}}-1)\sigma^{2}_{u}}$
(19)

where
$M=(\alpha_{b}^{2}+\sigma^{2}_{r})tr(\mathbf{H}_{RU}\mathbf{H}_{RU}^{h})^{-1}$.

####  IV-2 BS power optimization

When the relay power is assumed to be constant in equation (17a) the BS power
optimization sub-problem can be expressed as

$\displaystyle\max_{P_{b}}\frac{\frac{W}{2}\sum_{k=1}^{K}log_{2}(1+\gamma_{u,k})}{\zeta_{b}P_{b}+\zeta_{r}P_{r}+P_{sys}}$
$\displaystyle s.t.\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ 0\leq P_{b}\leq P_{b_{max}},$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ R_{k}\geq
R_{min,k},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ k=1,...,K.$ (20)

This sub-problem is similar to the problem of optimal relay power allocation.
The EE function in this case is also quasi-concave and therefore the PBA for
BS power variable can be used to find the optimal solution. In this case,
Algorithm 1 can be used by replacing index $b$ by index $r$ and vice-versa.
Using the same approach of previous part, $\rm{P_{b_{min}}}$ can be expressed
as

$P_{b_{min}}=\max_{k=1,...,K}\frac{(2^{2R_{min,k}}-1)M(\alpha^{2}_{r}+\sigma_{r}^{2}+\alpha^{2}_{u})}{\alpha^{2}_{r}}.$
(21)

1:$\rm{P_{r_{min}},P_{b_{min}},\leavevmode\nobreak\
P_{r_{max}},P_{b_{max}},(\epsilon_{b},\epsilon_{r})>0,\kappa}$

2:$\rm{\eta_{EE}^{*},\leavevmode\nobreak\ P_{r}^{*},\leavevmode\nobreak\
P_{b}^{*}}$

3:$\rm{initialize\leavevmode\nobreak\ P_{r}^{(0)},\leavevmode\nobreak\
P_{b}^{(0)},\leavevmode\nobreak\ l=0}$

4:$\rm{calculate\leavevmode\nobreak\ \eta_{EE}^{(l)}\leavevmode\nobreak\
from\leavevmode\nobreak\ \eqref{eq_ee}\leavevmode\nobreak\
for\leavevmode\nobreak\ P_{r}^{0}\leavevmode\nobreak\ ,P_{b}^{0}}$

5:repeat

6: $\rm{calculate\leavevmode\nobreak\ P_{r}^{*}\leavevmode\nobreak\
for\leavevmode\nobreak\ fixed\leavevmode\nobreak\
P_{b}^{(l)}\leavevmode\nobreak\ from\leavevmode\nobreak\
\eqref{eq_maxEE_SubRelay_Dl}\rightarrow\leavevmode\nobreak\ P_{r}^{(l+1)}}$

7: $\rm{update\leavevmode\nobreak\ P_{b_{min}}^{(l+1)}}$ as decribed in
Section IV Part 1

8: $\rm{calculate\leavevmode\nobreak\ P_{b}\leavevmode\nobreak\
for\leavevmode\nobreak\ fixed\leavevmode\nobreak\
P_{r}^{(l+1)}\leavevmode\nobreak\ from\leavevmode\nobreak\
\eqref{eq_maxEE_SubBS_Dl}\rightarrow\leavevmode\nobreak\ P_{b}^{(l+1)}}$

9: $\rm{update\leavevmode\nobreak\ P_{r_{min}}^{(l+1)}}$ as decribed in
Section IV Part 2

10: $\rm{update\leavevmode\nobreak\ l=l+1}$

11: $\rm{calculate\leavevmode\nobreak\ \eta_{EE}^{(l)}\leavevmode\nobreak\
from\leavevmode\nobreak\ \eqref{eq_ee}\leavevmode\nobreak\
for\leavevmode\nobreak\ P_{r}^{l}\leavevmode\nobreak\ ,P_{b}^{l}}$

12:until
$\rm{|\eta_{EE}^{(l)}-\eta_{EE}^{(l-1)}|\leq\epsilon\leavevmode\nobreak\
\leavevmode\nobreak\ }$

13:$\rm{P_{r}^{*}=P_{r}^{(l)},P_{b}^{*}=P_{b}^{(l)}\leavevmode\nobreak\
and\leavevmode\nobreak\ \eta_{EE}^{*}=\eta_{EE}^{(l)}}$

Algorithm 2 AOP Optimization for relay & BS Power

The PB-ODS algorithm can mitigate the complexity of the solution compared to
exhaustive search in both dimensions but the algorithm can still be considered
computationally complex since the search for one of optimal power values is
the exhaustive search. For further complexity reduction, it is proposed to use
the PBA algorithm alternatively for both dimensions until the convergence of
the algorithm which is denoted by AOP algorithm as described in algorithm 2.
As it can be seen, at each iteration the optimum value for one of the power
variables is computed using the PBA while keeping the other power variable
fixed and then this process is repeated for the other power variable.
According to this algorithm, the relay and BS power in the first part of the
algorithm are initialized and the corresponding EE is calculated. In each sub
problem, the problem of optimal power is solved individually for each power
variable according to the Algorithm 1 and after the optimization of both power
variables, the convergence condition is checked as the stopping criteria of
the algorithm. Finally, when the EE value converges, the BS and relay
transmission powers are determined from the sub-optimal value of their
corresponding problem and the near optimal EE value of the system is also
computed from the power variables.

##  V Simulation Results

In this section the simulation results are presented and discussed. The
simulation parameters of the system model are listed at table I. According to
the table, users are uniformly distributed at the edge of a circular cell
(800-900 meters from the center of the cell) in a 60-degree sector. Accoring
to the Long Term Evolution (LTE) standard , the bandwidth of each sub-channel
is set to be 180 KHz. The QoS constraint of all users are considered equal to
$1bps/Hz$.

Table I: Simulation parameters Parameter |  Value |  Parameter |  Value
---|---|---|---
Cell radius |  $900$ m |  $K$ |  $10$
User distance |  $800-900$ m |  Relay distance |  $425$ m
$\eta$ |  $3.5$ |  $P_{0}$ |  $27$ dBm
$P_{tx}$ |  $24$ dBm |  $P_{cod}/P_{dec}$ |  $29$ dBm
$P_{rx}$ |  $19$ dBm |  $T_{c}$ |  $32$ msec
$N_{0}$ |  $-100$ dBm |  $W$ |  $180$ KHz
$P_{b_{max}}$ |  $29$ dBm |  $R_{min}$ |  $1$ bps/Hz
$P_{r_{max}}$ |  $23$ dBm |  $\alpha$ |  $2$
$P_{syn}$ |  $28$ dBm |  $\zeta_{b},\zeta_{r}$ |  $0.3$
$L$ |  $10^{9}$ op/J |  $\kappa$ |  $1.1$

The following optimization methods are investigated.

  * •

Exhaustive search in two dimensions (TDS)

  * •

PBA for relay power and exhaustive search for BS power (PBR-ODSB)

  * •

PBA for BS power and exhaustive search for relay power (PBB-ODSR)

  * •

AOP in two dimensions (AOP)

Figure 2: Comparison of different algorithms on energy efficiency.

Fig. 2 shows the EE of the proposed system model for different power
allocation algorithms. Also, the performance of the PB-ODS algorithm
approaches the performance of the exhaustive search method. It can be
concluded that the use of the PBA algorithm in one dimension will be suitable
for finding the optimal value of the power. PB-ODS algorithm still uses an
exhaustive search method in one of two dimensions that has a high complexity.
Although the AOP method achieves lower performance than other algorithms, it
is extremely less complicated than other algorithms, which makes it suitable
for the cases that the processing speed of the algorithms is important. It can
also be seen that the performance of PBR-ODSB is better that PBB-ODSR. In
other words if the system has the process ability to use exhaustive search in
one of the dimensions, it is better to use exhaustive search for relay power
and PBA for BS power.

Fig. 3 represents the SE of the system in different modes as they are
described in table II. It can be seen that in the case of using maximum
transmission power for BS or relay, more SE is achieved for the system.
However, the aim of the optimization problem is to search for the BS and relay
transmission powers in order to maximize the EE of the system while the SE is
hold above a specified threshold.

Table II: Different power allocation modes Description | Name
---|---
$P_{b}=P_{b_{max}}\leavevmode\nobreak\ ,\leavevmode\nobreak\ P_{r}$ from PBA | Max BS
$P_{b}=P_{b_{min}}\leavevmode\nobreak\ ,\leavevmode\nobreak\ P_{r}$ from PBA | Min BS
$P_{r}=P_{r_{max}}\leavevmode\nobreak\ ,\leavevmode\nobreak\ P_{b}$ from PBA | Max relay
$P_{r}=P_{r_{min}}\leavevmode\nobreak\ ,\leavevmode\nobreak\ P_{b}$ from PBA | Min relay
$P_{r}\leavevmode\nobreak\ \&\leavevmode\nobreak\ P_{b}$ from PBA | AOP
$P_{r}\leavevmode\nobreak\ \&\leavevmode\nobreak\ P_{b}$ from Full Search | TDS
Figure 3: The spectral efficiency of the downlink transmission for different
power allocation modes (maximum power, minimum power, optimal power and sub-
optimal power).

Fig. 4 shows the EE of the system for different power allocation modes of
table II. Due to the usage of ZF pre-coding in the relay and the BS, there is
no interference among users. Therefore, it can be seen in the Fig. 3 that
increasing the transmission power of relay results in increasing the SE.
However, it can be seen that the EE of the system with maximum power
allocation of the relay and BS is less than optimal and sub-optimal
approaches. It should be again noted that the aim of this research is to
maximize the EE value which may result in SE performance degradation of the
system.

Figure 4: Energy efficiency of downlink transmission for different power
allocation modes.

##  VI Conclusion

In this paper the EE of a massive MIMO single-cell relay system in downlink
transmission is investigated. The relay is also assumed to be a massive MIMO
antenna which is a less investigated assumption in the literature. The goal is
to serve cell-edge users that face high pathloss. The problem is to maximize
the EE through optimal power allocation of relay and BS. The optimal solution
can be achieved using the exhaustive search algorithm which has a very high
computational complexity. Considering the quasi-concave behavior of the EE
function with respect to each one of power variables, we proposed to use PBA
in order to find the optimal value of the EE function with respect to each
power variable. Therefor, we suggest to use the PBA for one of power variables
following the exhaustive search method for the other power variable which is
called the PB-ODS algorithm. The simulation results show that the performance
of the PB-ODS algorithm approaches the performance of exhaustive search
algorithm with less amount of computational complexity. In order to further
reduce the complexity, it is suggested to use the PBA for both power variables
alternatively in an iterative manner which is called the AOP algorithm.
Simulation results shows that performance of the AOP algorithm is slightly
lower than the high complex optimal exhaustive algorithm but with a dramatic
decrease in computational complexity which makes it suitable for practical
scenarios. Also it is shown that due to the quasi concavity behavior of the EE
function, maximum power transmission does not necessarily results in the
maximum EE of the system.

## References

  * [1] A. Damnjanovic, J. Montojo, Y. Wei, T. Ji, T. Luo, M. Vajapeyam, T. Yoo, O. Song, and D. Malladi, “A survey on 3gpp heterogeneous networks,” _IEEE Wireless communications_ , vol. 18, no. 3, 2011.
  * [2] M. H. Bahonar and M. J. Omidi, “Centralized qos-aware resource allocation for d2d communications with multiple d2d pairs in one resource block,” in _Electrical Engineering (ICEE), Iranian Conference on_. IEEE, 2018, pp. 643–648.
  * [3] H. Gao, T. Lv, X. Su, H. Yang, and J. M. Cioffi, “Energy-efficient resource allocation for massive mimo amplify-and-forward relay systems.” _IEEE Access_ , vol. 4, pp. 2771–2787, 2016.
  * [4] B. Wang, J. Zhang, and A. Host-Madsen, “On the capacity of mimo relay channels,” _IEEE Transactions on Information theory_ , vol. 51, no. 1, pp. 29–43, 2005.
  * [5] H. Bolcskei, R. U. Nabar, O. Oyman, and A. J. Paulraj, “Capacity scaling laws in mimo relay networks,” _IEEE Transactions on Wireless Communications_ , vol. 5, no. 6, pp. 1433–1444, 2006.
  * [6] H. Cui, L. Song, and B. Jiao, “Multi-pair two-way amplify-and-forward relaying with very large number of relay antennas.” _IEEE Transactions Wireless Communications_ , vol. 13, no. 5, pp. 2636–2645, 2014.
  * [7] S. Jin, X. Liang, K.-K. Wong, X. Gao, and Q. Zhu, “Ergodic rate analysis for multipair massive mimo two-way relay networks,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 3, pp. 1480–1491, 2015.
  * [8] M. M. Naghsh, M. Soltanalian, P. Stoica, M. Masjedi, and B. E. Ottersten, “Efficient sum-rate maximization for medium-scale mimo af-relay networks,” _IEEE Trans. Wireless Communications_ , vol. 15, no. 9, pp. 6400–6411, 2016\.
  * [9] X. Jia, P. Deng, L. Yang, and H. Zhu, “Spectrum and energy efficiencies for multiuser pairs massive mimo systems with full-duplex amplify-and-forward relay.” _IEEE Access_ , vol. 3, pp. 1907–1918, 2015.
  * [10] D. P. Kudathanthirige and G. A. A. Baduge, “Multicell multiway massive mimo relay networks,” _IEEE Transactions on Vehicular Technology_ , vol. 66, no. 8, pp. 6831–6848, 2017.
  * [11] Z. Fang, X. Yuan, and X. Wang, “Towards the asymptotic sum capacity of the mimo cellular two-way relay channel,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 16, pp. 4039–4051, 2014.
  * [12] D. Kudathanthirige and G. A. A. Baduge, “Massive mimo configurations for multi-cell multi-user relay networks,” _IEEE Transactions Wireless Communications_ , vol. 17, no. 3, pp. 1849–1868, 2018.
  * [13] S. Huang, H. Chen, J. Cai, and F. Zhao, “Energy efficiency and spectral-efficiency tradeoff in amplify-and-forward relay networks,” _IEEE Transactions on Vehicular Technology_ , vol. 62, no. 9, pp. 4366–4378, 2013.
  * [14] J. Chen, H. Chen, H. Zhang, and F. Zhao, “Spectral-energy efficiency tradeoff in relay-aided massive mimo cellular networks with pilot contamination,” _IEEE Access_ , vol. 4, pp. 5234–5242, 2016.
  * [15] Z. Fang, W. Ni, F. Liang, P. Shao, and Y. Wu, “Massive mimo for full-duplex cellular two-way relay network: A spectral efficiency study,” _IEEE Access_ , vol. 5, pp. 23 288–23 298, 2017.
  * [16] J. Hoydis, S. Ten Brink, and M. Debbah, “Massive mimo in the ul/dl of cellular networks: How many antennas do we need?” _IEEE Journal on selected Areas in Communications_ , vol. 31, no. 2, pp. 160–171, 2013.