# A Multi-state Markov Model to Infer the Latent Deterioration Process From
the Maintenance Effect on Reliability Engineering of Ships

Hyunji Moon
IEOR Department
Columbia University, New York
<EMAIL_ADDRESS>
&Jungin Choi
Department of Statistics
Seoul National University
<EMAIL_ADDRESS>
&Seoyeon Cha
Department of Plant Science
Seoul National University
<EMAIL_ADDRESS>

###### Abstract

Maintenance optimization of naval ship equipment is crucial in terms of
national defense. However, the mixed effect of the maintenance and the pure
deterioration processes in the observed data hinders an exact comparison
between candidate maintenance policies. That is, the observed data-annual
failure counts of naval ships reflect counteracting actions between the
maintenance and deterioration. The inference of the latent deteriorating
process is needed in advance for choosing an optimal maintenance policy to be
carried out. This study proposes a new framework for the separation of the
true deterioration effect by predicting it from the current maintenance effect
through the multi-state Markov model. Using an annual engine failure count of
99 ships in the Korean navy, we construct the framework consisting of
imputation, transition matrix design, optimization, and validation. The
hierarchical Gaussian process model is used for the imputation and the three-
state Markov model is applied for the estimation of parameters in the
deterioration and maintenance effect. To consider the natural (deterioration)
and artificial (maintenance) effect respectively, the Bayesian HMM model with
a categorical distribution is employed. Computational experiments under
multiple settings showed the robustness of the estimated parameters, as well
as an accurate recovery of the observed data, thereby confirming the
credibility of our model. The framework could further be employed to establish
a reliable maintenance system and to reduce an overall maintenance cost.

_K_ eywords Multi-state Markov model $\cdot$ equipment reliability and
maintenance $\cdot$ optimization $\cdot$ imputation

##  1 Introduction

A maintenance policy is crucial both in terms of the safety and efficiency in
managing naval ships’ equipment. Maintenance policy includes several
controllable variables to be determined, such as an inspection frequency or an
acceptable maintenance standard. Too long inspection interval or overly
lenient standard for the repair would result in an unstable system, and the
resulting failure costs will be markedly increased. On the other hand, strict
maintenance with frequent inspections and excessively conservative standards
would yield an excellent budget waste.

For effective policy development, precise diagnosis of the current system
should be first addressed based on the historical failure data. For example,
annual engine failure count data is a significant source of information on two
main processes: deterioration and maintenance. However, the failure count data
ostensibly shows the combined effect; thus, identifying a hidden deteriorating
process is required.

Motivated by the fact that engine deterioration is a gradual process and its
future depends on the current state, we apply the multi-state Markov model to
predict deterioration. The multi-state model is defined as a stochastic
process that experiences a few possible states [5]. This stochastic model
could reflect much uncertainty, and the future state could be predicted with
randomness. Also, the multi-state approach can provide more interpretability
and an effective management policy in the real world since it perceives the
whole system as a gradual process between the perfect state (normal) and
complete failure (failure) [4, 6]. Moreover, computational complexity is
significantly reduced due to its discrete state modeling. Markovian multi-
state models have been widely used in infrastructure management for the
prediction of future deterioration states, such as bridge [7], healthcare
systems [8], and pavement [10] systems. However, the existing literature
mainly focused on risk assessment for a marine domain, and few have been on
the failure prediction [9, 1]. Many time-based multi-state Markov approaches
assume that the lifetime distribution is known, which generally does not hold
[11]. Since a regular investigation and maintenance have been implemented on
each system, most available failure data result from the current management
policy. Therefore, deterioration has been ignored in many past studies even
though the true deterioration process cannot be seen directly from observed
data. Reliable construction of transition probability matrix requires
consecutive deterioration data without any maintenance intervention [12]. To
deal with this problem, we predict the effect of the latent deterioration
process from the current maintenance policy.

This paper aims to suggest an overall framework for reliability engineering on
ships. At the same time, the data is incomplete, and the underlying
deterioration process is not directly observed due to the effect of the
current maintenance policy. Using the collected annual engine failure count
data for the recent ten years from the Korean navy, we estimate the true
deterioration process out of maintenance intervention by taking a nonlinear
optimization approach. Naval ship equipment failures are assumed as stochastic
degradation processes, and they are modeled as three multi-state continuous-
time Markov models. The transition probabilities between each deterioration
state are represented as a 3×3 transition probability matrix (TPM), which is
to be estimated. [13] computed the TPM by constructing the deterioration
distribution; however, they lacked the consideration of an age factor. This
study estimated the age-dependent transition probabilities with the transition
rate matrix obtained by the Kolmogorov equation, an inhomogeneous Markov
model. The results have shown different deterioration patterns for each age
period.

Figure 1: Process of the proposed model.

Figure 1 shows the overall process of our proposed model. First, the missing
data is imputed with the Gaussian process model to aid the next step,
parameter estimation of transition and maintenance matrix. Without this step,
estimation was highly inconsistent. Second, based on the observed data, the
best parameter representation of the deterioration and maintenance process was
designed. Then these parameters are estimated for given train sets. Lastly,
with the resulting parameter values, state series are predicted and compared
with real observed states. The last step is necessary for showing the validity
of our model.

The outline of this paper is as follows. In Section 2, we introduce our
overall modeling framework and multi-state Markov model. We formally define
the deterioration and maintenance matrix in Section 3 and also describe our
optimization process to separate the pure effect of the deterioration process
and estimate the hidden deterioration rate. Section 4 shows the result of our
approach, estimated deterioration parameters, and prediction accuracy. Lastly,
limitations and future developments are discussed in Section 5.

##  2 Problem setting

Annual failure counts of 99 naval ships for the recent ten years (2010-2019)
were retrieved from the equipment maintenance information system of the Korean
navy. 99 naval ships are classified into five engine types, and Figure 2
presents an overview of our data. The data has 99 vectors with a length of 31,
which represents the number of ships and their lifetime. The cells are colored
with different colors according to their engine type, and the uncolored cells
represent missing data. We arranged the failure data of 99 propulsion ship
engines according to their engine types (1 to 5). As we can see, the amount of
data for each engine type is highly imbalanced. Considering that the data was
recorded in the recent 10 years, the recorded lifetime region differs
according to its introduction time. For new models introduced after 2010
(e.g., engine type 5), failure counts were only recorded on younger ages (0-9
yrs), whilst old model’s data introduced in 1985 were constrained on older
ages (25-31 yrs). Moreover, a similarity between the data under the same
categories could be inferred; for example, same engine types tend to share
comparable age ranges and a scale of failure counts. Due to the security
problem of military data, only its scaled version is reported throughout the
paper.

Figure 2: Overview of the failure counts from the 99 Korean naval ships.

##  3 Imputation using a Gaussian Process

The original failure counts of 99 naval ships are largely missing and unevenly
spaced, imputation is required to achieve a robust multi-state Markov model.
After comparing several models including a hierarchical spline model and a
Poisson process, a hierarchical Gaussian process model was chosen for the
imputation. The hierarchical structure of the engine types in the system
validates the effectiveness of the applied hierarchical model [3].

We assume that failure counts follow a normal distribution with mean
$\mu_{t,j}$ and variance $\sigma^{2}_{k,j}$ for $j$th ship and time $t$, where
$k[j]$ is engine index for $j$th ship. This is engine index for jth ship. This
is based on the observation that the variances of failure counts between each
engine type look quite distinct. For example, while engine types 4 and 5 show
a large variance, type 3 shows a small variance of the given data.

$y_{t,j}\sim\mathcal{N}(\mu_{t,j},\sigma^{2}_{k[j]})$ (1)

Mean $\mu_{t,j}$ is defined as an additive Gaussian process which is the sum
of variables with normal distribution and Gaussian processes. $\mu$ is an
overall mean of $\mu_{t,j}$, and
$\theta_{t}^{age},\theta_{j}^{ship},and\theta_{k[j]}^{engine}$ follow a normal
distribution with zero mean and variance
$\sigma^{2}_{age},\sigma^{2}_{ship},and\sigma^{2}_{engine}$, respectively.

$\begin{split}\mu_{t,j}&=\mu+\theta_{t}^{\ age}+\theta_{j}^{\
ship}+\theta_{k[j]}^{\ engine}+\gamma_{t,j}+\delta_{t,\ k[j]}\\\ t&=1,2\dots
T,\ \ j=1,2,\dots N\\\ \end{split}$ (2)

$\begin{split}\theta_{t}^{\ age}&\sim\mathcal{N}(0,\sigma^{2}_{age}\ )\\\
\theta_{j}^{\ ship}&\sim\mathcal{N}(0,\sigma^{2}_{ship}\ )\\\ \theta_{k}^{\
engine}&\sim\mathcal{N}(0,\sigma^{2}_{engine}\ )\end{split}$ (3)

$\begin{split}\gamma_{j}&\sim\mathcal{N}(\textbf{0},\boldsymbol{K}_{l^{\gamma},{\alpha}^{\gamma}})\\\
\delta_{k}&\sim\mathcal{N}(\textbf{0},\boldsymbol{K}_{l^{\delta},{\alpha}^{\delta}})\end{split}$
(4)

$\gamma_{j}$ and $\delta_{k}$ are Gaussian processes whose $t$th elements are
$\gamma_{t,j}$ and $\delta_{t,k[j]}$, respectively. The covariance kernel for
the Gaussian processes $\gamma_{j}$ and $\delta_{k}$ are defined as an
exponentiated quadratic kernel. The exponentiated quadratic kernel defines the
covariance of each Gaussian process between $f(x_{i})$ and $f(x_{j})$ where
$f:\mathbb{R}^{D}\rightarrow\mathbb{R}$ as a function of the squared Euclidian
distance between $x_{i}\in\mathbb{R}^{D}$ and $x_{j}\in\mathbb{R}^{D}$ :

$\begin{split}\text{Cov}(f(x_{i}),f(x_{j}))&=k(x_{i},x_{j})\\\
&=\alpha^{2}\exp\left(-\frac{1}{2l^{2}}\sum_{d=1}^{D}(x_{i,d}-x_{j,d})^{2}\right)\end{split}$
(5)

with $\alpha$ and $l$ constrained to be positive.

Since the covariance function $\boldsymbol{K}_{l^{\gamma},{\alpha}^{\gamma}}$
and $\boldsymbol{K}_{l^{\delta},{\alpha}^{\delta}}$ follow an exponentiated
quadratic kernel, $ij$th element of
$\boldsymbol{K}_{l^{\gamma},{\alpha}^{\gamma}}$ and
$\boldsymbol{K}_{l^{\delta},{\alpha}^{\delta}}$ are defined as 6

$\begin{split}(\boldsymbol{K}_{l^{\gamma},{\alpha}^{\gamma}})_{ij}&={{\alpha}^{\gamma}}^{2}\exp\left(-\frac{1}{2{l^{\gamma}}^{2}}\sum_{d=1}^{D}(x_{i},d-x_{j},d)^{2}\right)\\\
(\boldsymbol{K}_{l^{\delta},{\alpha}^{\delta}})_{ij}&={{\alpha}^{\delta}}^{2}\exp\left(-\frac{1}{2{l^{\delta}}^{2}}\sum_{d=1}^{D}(x_{i},d-x_{j},d)^{2}\right)\end{split}$
(6)

The prior distribution for $l^{\gamma}$ and $l^{\delta}$ is the Weibull
distribution with shape parameter $k$ and scale parameter $\lambda$. Since
their mean largely depends on the scale parameter, we reparametrized
$l^{\gamma}$ and $l^{\delta}$ so that the scale parameter $\lambda$ is fixed
to 1.

$\begin{split}l_{\gamma}&\sim\text{Weibull}(k_{\gamma},\lambda_{\gamma})\\\
l_{\delta}&\sim\text{Weibull}(k_{\delta},\lambda_{\delta})\end{split}$ (7)

are reparametrized into

$\begin{split}l_{\gamma}^{\ s}&\sim\text{Weibull}(k_{\gamma},1),\ \ \ \
l_{\gamma}=l_{\gamma}^{s}*\lambda_{\gamma}\\\ l_{\delta}^{\
s}&\sim\text{Weibull}(k_{\delta},1),\ \ \ \
l_{\delta}=l_{\delta}^{s}*\lambda_{\delta}\end{split}$ (8)

Using this model, we gained an appropriate posterior distribution for each
parameter of the Gaussian model, thus imputing the missing data using the best
results in terms of prediction accuracy.

##  4 Separating the Deterioration and Maintenencnce effect using CTMC

###  4.1 Continuous time multi-state Markov model

In this study, we propose a three-state continuous-time Markov model to
describe the deterioration states over the naval ship’s life cycle as Figure
3. The ship’s deterioration state could be classified into three states
according to their annual failure counts by using a quantile classification
method as shown in Table 1. The corresponding failure counts are standardized
by a usual procedure, subtracting the mean and dividing by its standard error.
Each deterioration state of engines has a range from 1 to 3, where 1
corresponds to a "normal" state, 2 corresponds to a "near failure" state, and
3 to a complete "failure" state of engines. The state of each ship at time t
is assigned by the standardized annual failure counts as shown in Table 1.

Figure 3: Deterioration rate between three states. Table 1: Deterioration state of ship equipment according to the standardized annual failure counts. State | Status Description | Annual failure counts
---|---|---
1 | Normal | [ -1.8304, -0.3340 )
2 | Near Failure | [ -0.3340, -0.0703 )
3 | Failure | [ -0.0703, 2.1273 ]

A continuous-time multi-state model is for a continuous-time stochastic
process where individuals can occupy a finite number of states. Let
$\left\\{Y(t)|t\in T\right\\}$ represent the deterioration state of the
process at time $t$. The transition between states follows a continuous-time
Markov process, and transition probabilities only depend on a present state.
For a time-homogeneous Markov chain, we can write the transition probability
function from state $i$ at time $s$ to state $j$ at time $t$ as
$p_{ij}(s,t)=Pr\left\\{Y(t)=j|Y(s)=i\right\\}$. The transition probabilities
between all possible pairs $(i,j)$ are represented by a $n\times n$ matrix
called the transition probability matrix $P(s,t)$, where $n$ is the number of
possible condition states as shown in 9.

$P(s,t)=\begin{pmatrix}p_{11}(s,t)&p_{12}(s,t)&p_{13}(s,t)\\\
p_{21}(s,t)&p_{22}(s,t)&p_{23}(s,t)\\\
p_{31}(s,t)&p_{32}(s,t)&p_{33}(s,t)\end{pmatrix}$ (9)

In a continuous-time Markov process, transition rate matrix $Q$ is also
introduced to explain instantaneous transition rates between states. The
transition rate matrix $Q(t)$, also known as an intensity matrix, has an
element $q_{ij}$ (for $i\neq j$) denoting the rate departing from $i$ and
arriving in state $j$. Diagonal elements $q_{ii}$ are defined such that
${\displaystyle q_{ii}=-\sum_{j\neq i}q_{ij}.}$, satisfying the row-wise sums
of the matrix to zero. In other words, the off-diagonal elements of $Q$
represent the rates governing the exponentially distributed variables that are
used to describe the amount of time that elapses before a particular type of
base substitution occurs.

Transition matrix $P(t)$ can be calculated using the Kolmogorov’s forward and
backward equations, $P^{\prime}(t)=P(t)Q$ and $P^{\prime}(t)=QP(t)$,
respectively. The solution is uniquely derived in terms of the transition rate
matrix $Q$ as 10 by the matrix exponential of $Qt$.

$\begin{split}P(t)&=\exp(Qt)\\\
&=\sum_{r=0}^{\infty}Q^{r}\frac{t^{r}}{r!}\end{split}$ (10)

###  4.2 Deterioration and maintenance matrix

To separate the effect of the deterioration and the maintenance, we first
define a deterioration matrix $D(t)$, which is a function at the age $t$.
Deterioration matrix $D(t)$ should be an upper-diagonal matrix reflecting the
reality that engines continuously undergo degradation through the life cycle.
Also, since the extent of the deterioration varies between different age
ranges, it is constructed in a time-inhomogeneous manner. Older engines are
prone to be readily weakened, and it’s reasonable to design a more severe
deterioration matrix $D(t)$ for them. Equation 11 shows the deterioration rate
matrix $Q(t)$ for the three-state Markov model. This transition rate matrix
requires three distinct parameters for each time $t$, and is also upper-
diagonal, considering that the deterioration process can only occur in a
forward direction.

$Q(t)=\begin{pmatrix}-(\lambda_{t,1}+\lambda_{t,2})&\lambda_{t,1}&\lambda_{t,2}\\\
0&-\lambda_{t,3}&\lambda_{t,3}\\\ 0&0&0\\\ \end{pmatrix}$ (11)

As in equation 12, transition probabilities of the deterioration matrix $D(t)$
can be calculated solely with the given rate matrix $Q(t)$, which comes from
the solution of Kolmogorov’s forward and backward equation. Each element in 3
× 3 deterioration matrix $D(t)$ is the deterioration probability, and
$D_{ij}(t)$ is the probability of transitioning from state $i$ to $j$.

$\begin{split}D(t)&=\exp(Qt)\\\ &=\sum_{r=0}^{\infty}Q^{r}\frac{t^{r}}{r!}\\\
&=I+\sum_{r=1}^{\infty}Q^{r}\frac{t^{r}}{r!}\end{split}$ (12)

The resulting exponential matrix can be expressed in a closed form, because
the deterioration matrix $D(t)$ is positive-definite, which makes possible for
$Q$ to be always factorized as $Q=UDU^{-1}$ using eigen-decomposition, where
$D$ is the diagonal matrix of eigenvalues and $U$ is the matrix whose columns
are the corresponding eigenvectors. Using the diagonal matrix $D$, the
deterioration matrix $D(t)$ for three-state Markov model is computed as [14].
$p_{ij}(t)$ is equal to the occurence probability from state $i$ to $j$.

$D(t)=\begin{pmatrix}p_{11}(t)=e^{-(\lambda_{1}+\lambda_{3})t}&p_{12}(t)=\frac{\lambda_{1}e^{-\lambda_{3}t}(1-e^{-(\lambda_{1}+\lambda_{2}-\lambda_{3})t})}{(\lambda_{1}+\lambda_{2}-\lambda_{3})}&1-p_{11}(t)-p_{12}(t)\\\
0&e^{-\lambda_{3}t}&1-e^{-\lambda_{3}t}\\\ 0&0&1\end{pmatrix}$ (13)

In the same way, we define a maintenance matrix $M$ as in 14, which is a
transition probability matrix with two parameters. The maintenance matrix $M$
is multiplied to the state probability vector whenever the maintenance is
performed. The construction of this maintenance matrix is based on several
circumstantial assumptions as follows. First, maintenance interval is once a
year. Second, an imperfect maintenance makes the transition probability from
state 2 to be divided into state 1 and 2, which is $p_{21}$ and $1-p_{21}$ for
certain probability. Likewise, the transition from state 3 could be
parameterized with $p_{31}$, $p_{32}$, and $1-p_{31}-p_{32}$. However, the
preliminary experiment results showed that $1-p_{31}-p_{32}$ term mostly
converged to 0, we assumed it as 0 and replaced the probability $p_{31}$,
$p_{32}$ with $p_{31}$ and $1-p_{31}$, discarding the redundant parameter
$p_{32}$.

$M=\begin{pmatrix}1&0&0\\\ p_{21}&1-p_{21}&0\\\
p_{31}&1-p_{31}&0\end{pmatrix}$ (14)

###  4.3 HMM for the estimation of deterioration matrix

The main objective of our paper is to infer the latent deterioration effect
from the observed data intervened by the maintenance effect. This can be
achieved by estimating the underlying deterioration probabilities in the
deterioration matrix $D(t)$, which could be parametrized by the deterioration
rate parameters. The estimation of deterioration rate parameters $\lambda$
remains to be the main challenge.

For the parameter estimation, we consider a hidden Markov model (HMM) with the
assumption that observed failure count $y_{t}$ at time $t$ is generated based
on the hidden deterioration state vector
$d(t)=(d_{t,1},d_{t,2},d_{t,3})^{\top}$, where $d_{t,j}$ is the probability of
state $j$ at time $t$, hence a 3-simplex vector. According to HMM, the output
of failure count $y_{t}$ follows a categorical distribution,
$y_{t}\sim\textit{categorical}(d(t))$, which reflects stochastic randomness of
nature.

With respect to the Markov property, the present deterioration state vector
$d(t)$ only depends on the past vector $d(t-1)$. Starting from the initial
deterioration state vector, $d(0)=(1,0,0)^{\top}$, the ships continuously
undergo the deterioration process and the maintenance simultaneously at each
time, which is modeled with the matrix $D(t)$ and $M$ as mentioned. Taken
together, the deterioration state vector at time $t$ is formulated as below:

$d(t)=D(t)\\{\prod_{i=1}^{t-1}MD(i)\\}d(0),\ \
y_{t}\sim\textit{categorical}(d(t))$

We took a Bayesian paradigm for the estimation of established parameters with
the probabilistic statistical programming language, Stan, for a flexible
Bayesian statistical analysis [15]. Stan takes advantage of the Hamiltonian
Monte Carlo (HMC) algorithm to implement the Markov chain Monte Carlo (MCMC)
sampling, which is very efficient and powerful to estimate the posterior
distribution and where we are struggling to attain.

In our problem, the target parameters are $\lambda_{1}$, $\lambda_{2}$,
$\lambda_{3}$, $p_{21}$, and $p_{31}$, hereafter denoted $\theta$ all at once,
and the likelihood function $p(y|\theta)$ is a categorical distribution. As no
specific priors are assigned in our code, Stan automatically used a uniform
prior to yielding a final posterior distribution $p(\theta|y)$. Stan code to
fit the aforementioned model is illustrated in detail below.

    data {
      int<lower=0> N; // numbmer of obs
      int<lower=0> T; // number of time bins = max age
      int<lower=1>S; // number of states
      int<lower=1> P; // number of periods for inhomogenous dtr.
      int<lower=0, upper =S> states[N]; //observed state
      int obs2time[N]; // map obs to time
      int<lower=0, upper=S> initial_state;
    }

    transformed data {
      vector[S] initial;
      vector<lower=0>[S] alpha;
      vector<lower=0>[S] beta;
      initial = rep_vector(0, S);
      initial[initial_state] = 1;
      alpha = rep_vector(3,3);
      beta = rep_vector(3,3);
    }

    parameters {
      real<lower=0> rate[S];
      real<lower=0, upper=1> p21;
      real<lower=0, upper=1> p31;
    }

    transformed parameters {
      matrix[S, S] Dtr; // Deterioration
      matrix[S, S] Mnt; // Maintenance
      simplex[S] latent_states[T];
      matrix[S, S] tmp_p;
      // Maintenance
      // is left stoch. matrix, transpose of eq.14 from the paper
      Mnt = [[1, p21, p31],
             [0, 1-p21, 1-p31],
             [0, 0, 0]];
      // Deterioration by period
      // is left stoch. matrix, transpose of eq.13 from the paper
      tmp_p[1,1] = exp(-rate[1]- rate[2]);
      tmp_p[2,1] = rate[1] * exp(-rate[3]) * (1-exp(-(rate[1]+ rate[2] - rate[3]))) /
                  (rate[1]+ rate[2] - rate[3]);
      tmp_p[3,1] = exp(-rate[3]);
      Dtr = [[tmp_p[1,1], 0, 0],
                [tmp_p[2,1], tmp_p[3,1], 0],
                [1 - tmp_p[1,1] - tmp_p[2,1], 1 - tmp_p[3,1], 1]];
      // Inhomogenous Dtr
      latent_states[1] = Dtr * initial;
      for (t in 2:T){
        latent_states[t] =  (Dtr * Mnt) *latent_states[t-1];
        //matrix_power((Dtr * Mnt), (t-1)) * initial;
      }
    }

    model {
      for (n in 1:N){
        states[n] ~ categorical(latent_states[obs2time[n]]);
      }
    }

##  5 Result

###  5.1 Estimated parameters

As we parameterized deterioration probabilities and the maintenance
probabilities, the deterioration matrix $D(t)$ is parameterized by the rate
matrix $Q$ including 3 rate parameters for each 4 period, totally 12
deterioration rate parameters. The maintenance matrix $M$ is parameterized by
2 parameters, $p$ and $q$. With the observed deterioration data, we split the
total data into training and test sets by 94:5 ratio for 1000 different
combinations.

Table 2 shows the estimated average values and standard deviation of rate
parameters $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ . On the whole,
$\lambda_{1}$ and $\lambda_{3}$ are quite larger than $\lambda_{2}$, implying
that transition from state 1 to 2, and 2 to 3 is more likely than from state 1
to 3. Figure 4 shows the probability density graph of the three rate
parameters: $\lambda_{1},\lambda_{2}$, and $\lambda_{3}$, displaying similar
shapes between each parameter with certain peak values.

Table 2: Estimated mean value and their standard deviation of rate parameters: $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$ $\lambda_{i}$ | Mean | SD
---|---|---
$\lambda_{1}$ | 0.679 | 0.161
$\lambda_{2}$ | 0.274 | 0.149
$\lambda_{3}$ | 0.649 | 0.287
Figure 4: Resulting probability density graph of rate parameters:
$\lambda_{1},\lambda_{2}$, and $\lambda_{3}$.

Figure 5 shows the estimated distribution of every component of the
deterioration matrix $D(t)$. Since $D_{(t)}$ expresses the transition
probabilities between 3 deterioration states, $D_{12}(t)$, $D_{13}(t)$ and
$D_{23}(t)$ are set to zero, and $D_{33}(t)$ are set to 1 for every period.
Each element of $D(t)$ has distinctive peaks on its probability density. Table
4 shows the estimated mean value and their standard deviation of components of
the deterioration matrix $D_{(t)}$. The deterioration probability from state 1
to 1,2, and 3 shows almost equal probability around 0.3. Conversely, the
transition probability from state 2 to 3 is almost 0.6, while the staying
probability on state 2 is almost half of it. This implies that the
deteriorating process is highly likely from state 2.

Figure 5: Estimated probability density of the deterioration matrix $D(t)$. Table 3: Estimated average value and standard deviation of deterioration matrix: $D(t)$. $D_{ij}$ | Mean | SD
---|---|---
$D_{11}$ | 0.387 | 0.0302
$D_{12}$ | 0 | 0
$D_{13}$ | 0 | 0
$D_{21}$ | 0.300 | 0.0294
$D_{22}$ | 0.545 | 0.1603
$D_{23}$ | 0 | 0
$D_{31}$ | 0.313 | 0.0280
$D_{32}$ | 0.455 | 0.1603
$D_{33}$ | 1 | 0

On the other hand, for the maintenance parameters, $p21$ and $p31$ showed an
average value of 0.787 and 0.794, respectively. Table 4 shows the estimated
average values and standard deviation of rate parameters $p_{21}$ and
$p_{31}$. This indicates the transition probability from state 2 to 1 and from
state 3 to 2 after the yearly maintenance, implying that engine’s states tend
to go to state 1 with high probability for both state 2 and 3 engines after
the yearly maintenance. Thus, we observed the fine quality of the current
maintenance and repair of the Korean navy. Similarly, figure 6 shows the
resulting sampling results of two maintenance parameters, $p_{21}$ and
$p_{31}$.

Table 4: Estimated mean value and their standard deviation of maintenance parameters: $p_{21}$ and $p_{31}$ $p_{ij}$ | Mean | SD
---|---|---
$p_{21}$ | 0.787 | 0.167
$p_{31}$ | 0.794 | 0.167
Figure 6: Resulting sampling results of the maintenance parameters: $p_{21}$
and $p_{31}$.

According to the estimated values of $p_{21}$ and $p_{31}$, Figure 7 shows the
resulting distribution of every component of the maintenance matrix $M$. Since
$M$ is composed of maintenance probabilities, $M_{12}$,$M_{13}$,$M_{23}$ and
$M_{33}$ are fixed to 0 and $M_{11}$ is fixed to 1.

Figure 7: Estimated probability density of the maintenance matrix $M$.

###  5.2 Predicted states

The main purpose of the model is to predict the deteriorating state of each
engine so that management could be calculated based on their predictions.
Figure 8 shows the observed ratio of the 99 naval ship’s deterioration states,
representing their observed probabilities with the size of the black circles.
Red circles indicate the predicted states for each time. Predictions fit well
with the observed state from three out of the sample ship’s engines in most
cases. For example, between years 10 and 20, the flat area which corresponds
to the flat area of a bathtub in previous studies, the observed states are
concentrated in the first two states and so are the predictions.

Figure 8: Predicted states and observed states over time. Size of point
reflects ratio of states over 99 ships.

One of the best ways to summarize the accuracy of the prediction is by
measuring the distance between the observed and predicted states. Equation 15
shows the error measure between two state series, where $N$ denotes the number
of predicted ships, $d_{i}(t)$ and $\hat{d}_{i}(t)$ is the observed and
predicted deterioration state of $i$th ship at time $t$, respectively. Note
that the value of N is 5 for the test set and 94 for the training set.

Figure 15 shows the histogram of MSE values between the observed and predicted
states for each training and test set. A total of 1000 tests have been
repeated and test sets are constructed by randomly selected five engines.
Though the training error (average 20.7) is smaller than the test error
(average 20.9), MSE values of both training and test sets were barely
different, thus confirming the validity of our multi-state Markov model with
HMM.

$MSE=\frac{1}{N}\frac{1}{31}\sum_{i=1}^{N}\sum_{t=1}^{31}(d_{i}(t)-\hat{d}_{i}(t))^{2}$
(15)

Figure 9: Histogram of the training and test sets MSE values.

###  5.3 Simulation Based Calibration: SBC

Since we estimated the rate parameters in the deterioration matrix $D(t)$, and
maintenance parameters in the maintenance matrix $M$, we should validate
whether the MCMC computation is correct or not. Since our model is complicated
in implementations and algorithms, we validated inferences from Bayesian
algorithms capable of generating posterior samples through Simulation-Based
Calibration: SBC[20].

Figure 10 shows the SBC histograms for each model parameter which are the
maintenance parameters $p_{21}$, $p_{31}$, and rate parameters $\lambda_{1}$,
$\lambda_{2}$, and $\lambda_{3}$. Since SBC histograms for $p_{31}$,
$\lambda_{1}$, and $\lambda_{2}$ are fairly uniformly distributed and their
distributions are mostly included in the blue zone, which identifies accurate
computation and consistencies in model implementations. However, SBC
histrograms for $p_{21}$ and $\lambda_{3}$ are not uniformly distributed that
infers inaccurate computation.

Figure 10: The SBC histogram for each model parameter: $p_{21}$, $p_{31}$,
$\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$.

Graphical summaries in Figure 10 also the nature of the problems of
computations in $p_{21}$ and $\lambda_{3}$[20].

Since the SBC histogram for $p_{21}$ is $\cap$-shaped, it indicates that the
computed posterior distribution is over-dispersed relative to the prior
distribution. This implies that on average, the computed $p_{21}$ posterior
will be wider than the true $p_{21}$ posterior. Since $p_{21}$ reflects
maintenance from State 2 to 1, the true maintenance effect will be narrower
than the computed posterior distribution.

Since the SBC histogram for $\lambda_{3}$ is right-skewed, it indicates that
the algorithm is strongly biased towards larger values of $\lambda_{3}$ in the
true posterior. This might be caused by the identifiability of $\lambda_{3}$
and the heterogeneous Markov chain. Since $\lambda_{3}$ reflects deterioration
from State 2 to 3, the algorithm is biased on stronger deterioration from
’Near Failure’ to ’Failure’ status in the true posterior.

##  6 Conclusion

The multi-state model has a strong potential in reliability analysis. We have
applied this framework to the Korean Naval ships’ engine failure data and
established a three-state Markov model on the engine deterioration process.
The difficulty in maintenance policy suggestion is the inference of the pure
deterioration matrix since the current implemented policy is reflected in the
given data. We separately estimated the deterioration and maintenance effect
using nonlinear optimization based on the imputed data using the hierarchical
Gaussian process. Our approach provides a basis for a reasonable policy
comparison. We expect maintenance policy suggestions based on our methodology
to be the main direction for further research.

The limitation or further development of this study is as follows. First, the
deterioration matrix D(t) could be learned hierarchically. Currently, D(t) of
each engine type (1 to 5) is learned in a no-pooling way but it could be
updated to partial pooling, in other words, in a hierarchical way. Moreover,
the rate between four separated time regions in each ship could be partially
pooled as well. Though parameters would be added from this additional
structure, as has been introduced in [2] and [3], a hierarchical model with
its advanced flexibility could increase the model’s accuracy. Second, prior
knowledge of the engines (deterioration and operation patterns) could be
reflected. Priors could greatly affect the performance of the model if
properly understood in the context of the entire Bayesian analysis, from
inference to prediction to model evaluation [16].

## References

  * [1] Wang, Siqi and Yin, Jingbo and Khan, Rafi Ullah. (2020). The Multi-State Maritime Transportation System Risk Assessment and Safety Analysis. Sustainability, 12(14):5728
  * [2] Richard McElreath.Statistical rethinking: A Bayesian course with examples in R and Stan. CRC press, 2020.
  * [3] Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin. Bayesian data analysis. CRC press, 2013.
  * [4] GQ Qiu, S Huang, LL Zhu, XH Su, and Y Chen. Risk assessment of multi-state bayesian network in an oilgathering and transferring system.Procedia Engineering, 130:1514–1523, 2015.
  * [5] Philip Hougaard. Multi-state models: a review.Lifetime data analysis, 5(3):239–264, 1999.
  * [6] Krzysztof Kolowrocki and Joanna Soszynska-Budny.Reliability and Safety of Complex Technical Systems and Processes: Modeling–Identification–Prediction-Optimization. Springer Science & Business Media, 2011.
  * [7] Maatern-Jan Kallen. Markov processes for maintenance optimization of civil infrastructure in the netherlands.2007.
  * [8] Jaime González-Domínguez, Gonzalo Sánchez-Barroso, and Justo García-Sanz-Calcedo. Scheduling of preventivemaintenance in healthcare buildings using markov chain.Applied Sciences, 10(15):5263, 2020.
  * [9] Shahrzad Faghih-Roohi, Min Xie, and Kien Ming Ng. Accident risk assessment in marine transportation viamarkov modelling and markov chain monte carlo simulation.Ocean engineering, 91:363–370, 2014.
  * [10] Abbas A Butt, Mohamed Y Shahin, Kieran J Feighan, and Samuel H Carpenter.Pavement performance predictionmodel using the Markov process. Number 1123. 1987.
  * [11] Bram de Jonge and Philip A. Scarf. A review on maintenance optimization.European Journal of Operational Research, 285(3):805 – 824, 2020.
  * [12] George Morcous. Performance prediction of bridge deck systems using markov chains.Journal of performance of Constructed Facilities, 20(2):146–155, 2006.
  * [13] Li Li, Lijun Sun, and Guobao Ning. Deterioration prediction of urban bridges on network level using markov-chain model.Mathematical Problems in Engineering, 2014, 2014.
  * [14] Edmund Jones, David Epstein, and Leticia García-Mochón. A procedure for deriving formulas to convert transition rates to probabilities for multistate markov models.Medical Decision Making, 37(7):779–789, 2017.
  * [15] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus ABrubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: a probabilistic programming language. Grantee Submission, 76(1):1–32, 2017.
  * [16] Andrew Gelman, Daniel Simpson, and Michael Betancourt. The prior can often only be understood in the contextof the likelihood.Entropy, 19(10):555, 2017.
  * [17] Killick, Rebecca and Beaulieu, Claudie and Taylor, Simon and Hullait, Harjit. Package ‘EnvCpt’. 2020.
  * [18] Yongfeng FANG, Wenliang TAO and Tee Kong Fah. Reliability analysis of multi-state engine units utilizing time-domain response data. Journal of Systems Science and Information, 4(4):354–364, 2016.
  * [19] Bent Natvig.Multistate Systems Reliability Theory with Applications. doi:10.1002/9780470977088, 2011.
  * [20] Sean Talts, Michael Betancourt, Daniel Simpson, Aki Vehtari, Andrew Gelman. Validating Bayesian Inference Algorithms with Simulation-Based Calibration. 2018.


# Study on Precoding Optimization Algorithms
in Massive MIMO System with Multi-Antenna Users

Evgeny Bobrova,b, Dmitry Kropotova,c, Sergey Troshinc and Danila Zaevb CONTACT
Evgeny Bobrov. Email<EMAIL_ADDRESS>aLomonosov MSU, Russia; bMRC, Huawei
Technologies, Russia; cHSE University, Russia

###### Abstract

The paper studies the multi-user precoding problem as a non-convex
optimization problem for wireless multiple input and multiple output (MIMO)
systems. In our work, we approximate the target Spectral Efficiency function
with a novel computationally simpler function. Then, we reduce the precoding
problem to an unconstrained optimization task using a special differential
projection method and solve it by the Quasi-Newton L-BFGS iterative procedure
to achieve gains in capacity. We are testing the proposed approach in several
scenarios generated using Quadriga — open-source software for generating
realistic radio channel impulse response. Our method shows monotonic
improvement over heuristic methods with reasonable computation time. The
proposed L-BFGS optimization scheme is novel in this area and shows a
significant advantage over the standard approaches. The proposed method has a
simple implementation and can be a good reference for other heuristic
algorithms in this field.

###### keywords:

Optimization, Massive MIMO, Precoding, SVD, L-BFGS, Interior-Point.

##  1 Introduction

Wireless channels with multiple input and multiple output (MIMO) provide
significantly more capacity than their counterparts with one input and one
output. Thus, the MIMO system is an essential technology in modern wireless
telecommunications, including Wi-Fi and 5G systems [2, 5]. It allows using
antenna array clusters to send multiple signal beams for multiple user devices
simultaneously. The proper construction of such beams is called beamforming or
precoding procedure [10, 14]. In the linear channel assumptions, the correct
precoding is a complex matrix of the general form with given constraints,
which corresponds to the physical limitations of the system. We measure the
quality of the obtained precoding using the well-known functions such as
Signal-to-Interference-and-Noise Ratio (SINR) [20] and Spectral Efficiency
(SE) [19].

The standard precoding algorithms, which are well-known in the literature, are
Maximum-Ratio Transmission (MRT) [12], Zero-Forcing (ZF) [24] and Regularized
ZF (RZF) [21] (including its recent variant — Adaptive RZF (ARZF) [3]). All
these algorithms have analytical formulas, without taking into account the
target function of Spectral Efficiency (SE) or do it implicitly, maximizing
the numerator of SINR using the MRT algorithm, or reducing the denominator of
SINR using the ZF algorithm. This leads to simple but non-optimal precoding
solutions. Therefore, one of our goal is to study the potential improvement of
precoding methods.

While during downlink (DL) on the transmitter side a base station apply
precoding matrices, symmetrically, on the receiver side users apply detection
matrices. The knowledge of the detection of the user allows constructing more
efficient precoding for the base station, some works even consider joint
construction of precoding and detection algorithms, see e.g. [17]. In general,
we do not have exact information about the user detection, but we can mean a
type of it. The most common detections are well-described in the literature,
such as Minimum MSE (MMSE) [13, 23] and Interference Rejection Combiner (MMSE-
IRC) [16], which can be applied in order to eliminate the multiuser
interference. We also utilize the recently proposed theoretical Conjugate
Detection (CD) [3].

In this article, we investigate the case when one transmitter (a base station)
sends data to several receivers (users), which is also called a downlink (DL)
procedure. We constrain ourselves to the case of linear transmitter and the
receiver. For the case when each user has only one receiving antenna under the
total power restriction at the base station, optimal linear precodings
involving Mimimum MSE (MMSE) [13, 23] function were derived in the papers [11,
1, 3]. In contrast, we consider the most practical formulation of the problem
with multiple-antenna users [18] and per-antenna power constraints at the base
station [22].

Thus, considered problem is constrained optimization programming task with
non-convex function and multiple constraints. It can be solved by using
suboptimal heuristic methods [12, 24], interior point methods [8], projection-
based methods [7], or any of them with reparametrization (variables
transform). In this work, we find a convenient reparametrization, that reduces
the constrained optimization problem to an unconstrained one. This reduction
allows to use an advanced unconstrained gradient-based methods. The proper
reparametrization is one of the main contributions of the article. Even with
the right reparameterization, the target function is computationally complex
and must be simplified to accelerate convergence. Proper simplification is the
second and new approach presented in this article.

In this study we, firstly, approximate the SE [19] objective function using
the Conjugate Detection CD [3] to simplify the optimization procedure as much
as possible. As the second step, we reduce the precoding problem to an
unconstrained optimization task and solve using the Quasi-Newton L-BFGS
iterative procedure [26] to achieve the maximum possible transmitting quality.
We write an end-to-end differentiable projection-based method, where gradients
are taken both for the functional and the projection, which tends to very fast
convergence of the proposed method. The problem to be solved is not convex,
and so we find the local maximum of the SE.

The novelty of the proposed method consists in simplifying the objective
function SE, which makes the task more attractive from a computational point
of view. Secondly, the reparametrization of precoding for an unconstrained
task was invented. We called this parametrization Differentiable projection.
In addition, classical parameterization using the Softmax function is also
described.

Our parametrization method is a know-how and is designed to solve the problem
of precoding. In fact, it can be generalized to other topics of constrained
optimization problems. In combination with the simplified function, we
proposed the target problem of unconstrained optimization as a reference
algorithm for evaluating the upper bound of quality. Finally, all our research
is based on realistic power per antenna power constrains and multiple antenna
users, which in itself is a separate problem for finding the right approach to
solving.

All investigated algorithms were tested in several scenarios of Quadriga [9] —
an open-source software for generating realistic radio channel impulse
response. Our method shows monotonic improvement over heuristic methods with
reasonable computational time.

We adopt commonly used notational conventions throughout the paper. Matrices
and vectors will be denoted by bold-face upper and lower case letters,
respectively. Furthermore, $(\cdot)^{\mathrm{H}}$ denotes Hermitian transpose.

The paper is organized as follows. In Section 2 we describe the system model
of the studied Massive MIMO network, introduce quality measures, types of
detection matrices, approximated quality function, power constraints and
problem statement. In Section 3 we study solutions for precoding problem,
reference methods, proposed optimization solution and its relative
computational complexity. In Section 4 we provide numerical results in the
Quadriga radio simulator. Section 5 contains the conclusion.

##  2 System Model

In 5G cellular networks, we look at a precoding problem in multi-user massive
MIMO communication. A single base station with multiple transmit antennas is
used. These antennas send signals to multiple users at the same time. Each
user also has multiple antennas for receiving signals. The transmission
quality between the base station and the users is measured by the base
station. The precoding problem entails determining an appropriate weighting
(phase and gain) for transmitting the signal in order to maximise the signal
power at the receivers’ output.

Figure 1: An example of the MIMO precoding usage. The problem is to find an
optimal precoding matrix $\bm W$ of the system given the target SE function
(3) and the per-antenna power constraints (8). Figure 2: An example of the
MIMO transmission system in the matrix form. Multi-User precoding $\bm W$
allows transmitting different information to various users simultaneously.

In our system we consider $K$ users and for $k$-th user we would like to
transmit him $L_{k}$ symbols. In total, we would like to transmit a vector
$\bm x\in\mathbb{C}^{L}$, where $L=L_{1}+\dots+L_{K}$. We multiply the
transmitted vector by a precoding matrix $\bm W\in\mathbb{C}^{T\times L}$,
where $T$ is a total number of transmit antennas on a base station. Then we
transmit the precoded symbols for the all users. Suppose that $k$-th user has
$R_{k}$ receive antennas and $\bm H_{k}\in\mathbb{C}^{R_{k}\times T}$ is a
channel between $k$-th user and the base station. Then $k$-th user receives
$\bm H_{k}\bm W\bm x+\bm n_{k}$, where $\bm n_{k}$ is a Gaussian noise. The
noise appears as a result of thermal distortions of the system. Finally,
$k$-th user applies a detection of transmitted symbols by multiplying the
received vector by a detection matrix $\bm G_{k}\in\mathbb{C}^{L_{k}\times
R_{k}}$. The final vector of detected symbols of the all users is denoted by
$\bm r\in C^{L}$. The whole process of transmitting symbols is presented in
Fig. 2 [3]. Note that the linear precoding and detection can be implemented by
simple matrix multiplications. Usually the total number of base station
antennas, user antennas and transmitted symbols are related as $L_{k}\leq
R_{k}\leq T$.

Finally, the Multi-User MIMO model is described using the following linear
system:

$\bm r=\bm G(\bm H\bm W\bm x+\bm n).$

Where $\bm r\in\mathbb{C}^{L}$ is a received vector, and $\bm
x\in\mathbb{C}^{L}$ is a sent vector, and $\bm H\in\mathbb{C}^{R\times T}$ is
a channel matrix, and $\bm W\in\mathbb{C}^{T\times L}$ is a precoding matrix,
and $\bm G\in\mathbb{C}^{L\times R}$ is a block-diagonal detection matrix, and
$\bm n\sim\mathcal{CN}(0,I_{L})$ is a noise-vector. The constant $T$ is the
number of transmit antennas, $R$ is the total number of receive antennas, and
$L$ is the total number of transmitted symbols in the system. They are
typically related as $L\leq R\leq T$. Each of the matrices $\bm G,\bm H,\bm W$
decomposes by $K$ users. The numbers $K,L,R,T$ and channel matrix $\bm
H\in\mathbb{C}^{R\times T}$ are predefined in the system.

It is convenient [18] to represent the user channel matrix $k$ via its reduced
Singular Value Decomposition (SVD):

$\bm H_{k}=\bm U^{\mathrm{H}}_{k}\bm S_{k}\bm V_{k},\quad\bm U_{k}\bm
U_{k}^{\mathrm{H}}=\bm U_{k}^{\mathrm{H}}\bm U_{k}=\bm I_{R_{k}},\quad\bm
S_{k}={\rm diag}\\{s_{1},\dots,s_{R_{k}}\\},\quad\bm V_{k}\bm
V^{\mathrm{H}}_{k}=\bm I_{R_{k}}.$

Where the channel matrix for user $k$, $\bm H_{k}\in\mathbb{C}^{R_{k}\times
T}$ contains channel vectors $\bm h_{i}\in\mathbb{C}^{T}$ by rows, the
singular values $\bm S_{k}\in\mathbb{C}^{R_{k}\times R_{k}}$ are sorted by
descending, $\bm U_{k}\in\mathbb{C}^{R_{k}\times R_{k}}$ is a unitary matrix
of left singular vectors, and matrix $\bm V_{k}\in\mathbb{C}^{R_{k}\times T}$
consists of right singular vectors.

Collecting all users together, we can write the following decomposition:

######  Lemma 2.1.

[3] Denote $\bm H=[\bm H_{1},\dots,\bm H_{K}]\in\mathbb{C}^{R\times T}$ the
concatenation of individual channel rows $\bm H_{k}$. The following
representation holds:

$\bm H=\bm U^{\mathrm{H}}\bm S\bm V.$

where the $\bm H\in\mathbb{C}^{R\times T}$, and $\bm S=\text{diag}(\bm
S_{k})\in\mathbb{C}^{R\times R}$, and $\bm U=\text{bdiag}(\bm
U_{k})\in\mathbb{C}^{R\times R}$ is a block-diagonal unitary matrix, but $\bm
V=[\bm V_{1},\dots,\bm V_{K}]\in\mathbb{C}^{R\times T}$ has a general form.

###  2.1 Quality Measures

The quality functions are based not on the actual sending symbols
$x\in\mathbb{C}^{L}$, but on some distribution of them [2]. Thus, we get the
common functions for all the assumed symbols, which can be sent using the
specified precoding. The Signal-to-Interference-and-Noise (SINR) functional of
the $l=l(k)$-th symbol and user k, is defined as:

$\textrm{SINR}_{l}(\bm W,\bm H_{k},\bm g_{l},\sigma,P)=\dfrac{|\bm g_{l}\bm
H_{k}\bm w_{l}|^{2}}{\sum_{i\neq l}^{L}|\bm g_{l}\bm H_{k}\bm w_{i}|^{2}+\|\bm
g_{l}\|^{2}\frac{\sigma^{2}}{P}},\quad\forall l\in\mathcal{L}_{k},$ (1)

where $\bm w^{l}\in\mathbb{C}^{T}$ is the $l$-th column of the precoding
matrix.

To get the formula for the $k$-th user SINR, namely the effective SINR, that
shows the signal quality of the user, we average his $L_{k}$ per-symbol SINRs
(1) by the geometric mean:

$\textrm{SINR}_{k}^{eff}(\bm W,\bm H_{k},\bm
G_{k},\sigma,P)=\Big{(}{\prod\nolimits_{l\in\mathcal{L}_{k}}\textrm{SINR}_{l}(\bm
W,\bm H_{k},\bm g_{l},\sigma,P)}\Big{)}^{\frac{1}{L_{k}}},$ (2)

where the detection matrix $G_{k}\in\mathbb{C}^{L_{k}\times R_{k}}$ contains
vector of the $k$-th user symbols $g_{l}$.

To get the final functional of Spectral Efficiency, we apply the Shannon’s
formula over all effective user SINRs (2):

$\textrm{SE}(\bm W,\bm H,\bm
G,\sigma,P)=\sum_{k=1}^{K}L_{k}\log_{2}(1+\textrm{SINR}_{k}^{eff}(\bm W,\bm
H_{k},\bm G_{k},\sigma,P))\rightarrow\max\limits_{\bm W}$ (3)

The formula (3) shows the throughput capacity of the system as a whole. Now it
depends on the whole channel and detection matrices $\bm
H\in\mathbb{C}^{R\times T}$ and $\bm G\in\mathbb{C}^{L\times R}$. We will use
it as our final score function, compare algorithms and find precoding matrices
$\bm W\in\mathbb{C}^{T\times L}$, solving the problem of maximizing this
function.

Additionally, we consider the function of Single-User SINR (SUSINR):

$\textrm{SUSINR}(\widetilde{\bm
S},\sigma^{2},P)=\frac{P}{\sigma^{2}}\bigg{(}\prod_{k=1}^{K}\frac{1}{L_{k}}\bigg{(}\prod\nolimits_{l\in\mathcal{L}_{k}}s_{l}^{2}\bigg{)}^{\frac{1}{L_{k}}}\bigg{)}^{\frac{1}{K}}.$
(4)

The formula (4) reflects the quality of the channel for the specified user
without taking into account the others. It depends on the greatest $L_{k}$
singular values $\widetilde{\bm S}_{k}\in\mathbb{R}^{L_{k}\times L_{k}}$ of
the $k$-th user channel matrix $\bm H_{k}\in\mathbb{C}^{R_{k}\times T}$.
Finally, the matrix $\widetilde{\bm S}\in\mathbb{C}^{L\times L}$ contains all
singular values on the main diagonal, such as $\widetilde{\bm
S}=\text{diag}(\widetilde{\bm S}_{1}\dots\widetilde{\bm
S}_{K})\in\mathbb{C}^{L\times L}$. We will use this function as a universal
channel characteristic in our experiments.

###  2.2 Detection Matrices

Detection aims to reverse the channel with precoding to an identity matrix. We
study the most common forms of the detection matrix $\bm G$, which is called
Minimum MSE (MMSE) [13], Interference Rejection Combine (MMSE-IRC) [16], and
Conjugate Detection (Conjugate) (CD). In our work we will focus our attention
on the MMSE-IRC and CD methods.

Firstly, MMSE detection realizes the following rule:

$\bm G^{\textrm{MMSE}}_{k}(\bm A_{k})=\bm A_{k}^{\mathrm{H}}\left(\bm A_{k}\bm
A_{k}^{\mathrm{H}}+\lambda\bm I\right)^{-1},\quad\bm A_{k}=\bm H_{k}\bm
W_{k}.$

The scalar $\lambda=\frac{\sigma^{2}}{P}$ is the system noise-to-signal ratio.
The constant $P$ refers to the base station power and $\sigma^{2}$ refers to
the system noise. Thus, the method implicitly consider the noise equal for all
symbols: $\bm n=\frac{\sigma^{2}}{P}\bm 1\in\mathbb{C}^{L}$, and detection
$\bm G$ which tends to eliminate it. One may notice that the use of the method
requires solving a system of linear equations of the size of layers $L$.

Secondly, MMSE-IRC detection uses additional information about noise-
covariance matrix $\bm R_{uu}^{k}$ as follows:

$\bm G^{\textrm{IRC}}_{k}(\bm A_{k},\bm R_{uu}^{k})=\bm A_{k}^{\mathrm{H}}(\bm
A_{k}\bm A_{k}^{\mathrm{H}}+\bm R_{uu}^{k}+\lambda\bm I)^{-1},\quad\bm
A_{k}=\bm H_{k}\bm W_{k},$ (5)

where the covariance matrix has the form as:

$\bm R_{uu}^{k}=\bm H_{k}\left(\bm W\bm W^{\mathrm{H}}-\bm W_{k}\bm
W_{k}^{\mathrm{H}}\right)\bm H_{k}^{\mathrm{H}}=\bm H_{k}\left(\sum_{u=1,u\neq
k}^{K}\bm W_{u}\bm W_{u}^{\mathrm{H}}\right)\bm H_{k}^{\mathrm{H}}.$

######  Lemma 2.2.

Using the linear property, you can get the following representation of the
detection matrix MMSE-IRC:

$\bm G^{\textit{IRC}}_{k}(\bm H_{k},\bm R_{uu}^{k})=\left(\bm H_{k}\bm
W_{k}\right)^{\mathrm{H}}(\bm H_{k}\bm W_{k}(\bm H_{k}\bm
W_{k})^{\mathrm{H}}+\bm R_{uu}^{k}+\lambda\bm I)^{-1}=\\\ =\left(\bm H_{k}\bm
W_{k}\right)^{\mathrm{H}}(\bm H_{k}\bm W_{k}(\bm H_{k}\bm
W_{k})^{\mathrm{H}}+\bm H_{k}\left(\sum_{u=1,u\neq k}^{K}\bm W_{u}\bm
W_{u}^{\mathrm{H}}\right)\bm H_{k}^{\mathrm{H}}+\lambda\bm I)^{-1}=\\\
=\left(\bm H_{k}\bm W_{k}\right)^{\mathrm{H}}(\bm H_{k}\bm W(\bm H_{k}\bm
W)^{\mathrm{H}}+\lambda\bm I)^{-1}.$

Thirdly, Conjugate Detection can be written in the following form:

$\bm G^{C}:=\widetilde{\bm S}^{-1}\widetilde{\bm U}\in\mathbb{C}^{L\times
R},\quad\bm G^{C}_{k}=\widetilde{\bm S}_{k}^{-1}\widetilde{\bm
U}_{k}\in\mathbb{C}^{L_{k}\times R_{k}}.$ (6)

Where $\widetilde{\bm S}_{k}\in\mathbb{C}^{L_{k}\times L_{k}}$ contains the
first $L_{k}$ largest singular values, and $\widetilde{\bm
U}_{k}\in\mathbb{C}^{L_{k}\times R_{k}}$ contains the first $L_{k}$
corresponding singular vectors. The Conjugate and MMSE-IRC detection matrices
are very closely related to each other, assuming low noise and ZF or RZF
Precoding, we provide a proof of this fact in the work [4]. We can use the
Conjugate detection at the base station to configure the precoding matrix
assuming MMSE-IRC for user equipments.

###  2.3 Approximated Quality Function

Assuming the conjugated detection $\bm G^{C}$ (6) we have come to the
following approximated function of the SINR [3]:

$\textrm{SINR}_{l}^{C}(\bm W,\widetilde{\bm
v_{l}},s_{l},\sigma^{2},P)=\frac{|\widetilde{\bm v_{l}}\bm
w_{l}|^{2}}{\sum_{i\neq l}^{L}|\widetilde{\bm v_{l}}\bm
w_{i}|^{2}+s_{l}^{-2}\frac{\sigma^{2}}{P}}$

Spectral Efficiency function can be simplified in the following way, where
$\widetilde{\bm v_{l}}\in\mathbb{C}^{T}$ is the singular vector of the $l$-th
symbol, and $\bm s_{l}\in\mathbb{R}$ is the singular value of the $l$-th
symbol:

$\textrm{SE}^{C}(\bm W,\widetilde{\bm V},\widetilde{\bm
S},\sigma,P)=\sum_{l=1}^{L}\log_{2}(1+\textrm{SINR}_{l}^{C}(\bm
W,\widetilde{\bm v_{l}},s_{l},\sigma,P))=\\\
\sum_{l=1}^{L}\log_{2}\Big{(}\sum_{i=1}^{L}|\widetilde{\bm v_{l}}\bm
w_{i}|^{2}+s_{l}^{-2}\frac{\sigma^{2}}{P}\Big{)}-\sum_{l=1}^{L}\log_{2}\Big{(}\sum_{i\neq
l}^{L}|\widetilde{\bm v_{l}}\bm
w_{i}|^{2}+s_{l}^{-2}\frac{\sigma^{2}}{P}\Big{)}\rightarrow\max\limits_{\bm
W}.$ (7)

One may notice that we have completely moved away from user antennas of shapes
$R_{k}$ and $R$ and work only with user layers of shapes $L_{k}$ and $L$. The
formula now does not depend on channel matrix $\bm
H_{k}\in\mathbb{C}^{R_{k}\times T}$ but on the eigenvectors of the $l$-th
layer $\bm v_{l}\in\mathbb{C}^{T}$, which has length $T$ by the number of
antennas. The inversed squared singular values $s_{l}^{-2}\in\mathbb{R}$ scale
the noise power.

###  2.4 Power Constraints

We formulate the realistic per-antenna power constraints [25]. Since we have
$T$ equal transmitter antennas, the power limitation applied to each of them
is $\frac{P}{T}$ . The antenna power can be described in terms of the row
norms of the precoding matrix:

$\|\bm w^{m}\|^{2}\leq\frac{P}{T}\;\forall m=1\dots T.$ (8)

It is clear that per-antenna constraints satisfies the total power: $\|\bm
w^{1}\|^{2}+\dots+\|\bm w^{T}\|^{2}\leq P$, which is the sum-power constraint
across all transmit antennas. While analytically attractive, such a sum-power
constraint is often unrealistic in practice. In a physical implementation of a
multi-antenna base station, each antenna has its own power amplifier in its
analog front-end and is limited individually by the linearity of the power
amplifier. Thus, a power constraint imposed on a per-antenna basis is more
realistic.

###  2.5 Problem Statement

We formulate the constrained smooth optimization problem: —l— W ∈C ^T ×LSE(W,
H, G, σ, P) ∥w^m∥^2≤PT, m = 1 …T Where $\bm w^{m}\in\mathbb{C}^{L}$ in the
$m$-th row of the precoding matrix.

##  3 Solutions for Precoding

###  3.1 Reference Methods

The precoding matrix $\bm W$ is responsible for the beamforming from the base
station to the users. The linear methods for precoding do the following.
Firstly, the linear solutions obtain singular value decomposition for each
user $\bm H_{k}=\bm U_{k}^{\mathrm{H}}\bm S_{k}\bm
V_{k}\in\mathbb{C}^{R_{k}\times T}$ (Lemma 2.1) and take the first $L_{k}$
singular vectors $\widetilde{\bm V}_{k}\in\mathbb{C}^{L_{k}\times T}$ which
attend to the first $L_{k}$ greatest singular values [18]. All these matrices
are concatenated to the one matrix $\widetilde{\bm V}\in\mathbb{C}^{L\times
T}$.

Finally, the precoding matrix is constructed from the obtained singular
vectors. We describe linear methods for constructing a precoding matrix. We
meet the power constraints using the scalar post-adjustment. We divide the
precoding matrix on its maximal row norm as $\max\\{\|\bm w^{t}\|\\}_{t=1\dots
T}$ and scale them on $\sqrt{\frac{P}{T}}$. Thus, the per-antenna power
constraints can satisfy the following condition: $\|\bm
w^{m}\|^{2}\leq\frac{P}{T}\;\forall m$.

Maximum-Ratio Transmission is the simplest way of computing the precoding
matrix. It takes the Hermitian conjugate of the channel singular vectors [12]:

$\bm W_{MRT}=\mu\widetilde{\bm V}^{\mathrm{H}}\bm P\in\mathbb{C}^{T\times L}$
(9)

Where matrix $\bm P\in\mathbb{R}^{L\times L}$ denotes power allocation between
symbols and constant $\mu>0$ utilizes the power constraints. In such a way,
the signal beams can be sent directly to users without considering their
interaction. The Maximum-Ratio approach is preferred in noisy systems where
the noise power is higher than inter-user interference.

Zero-Forcing is the next modification of the precoding algorithm. which
performs decorrelation of the symbols using inverse correlation matrix of the
channel vectors [24]:

$\bm W_{ZF}=\mu\widetilde{\bm V}^{\mathrm{H}}(\widetilde{\bm V}\widetilde{\bm
V}^{\mathrm{H}})^{-1}\bm P\in\mathbb{C}^{T\times L}.$ (10)

With matrix $\bm P\in\mathbb{R}^{L\times L}$ and constant $\mu>0$, which have
the same meaning as in the previous section. Such precoding construction sends
the signal beams to the users without creating any interference between them.
Different from the previous method, the Zero-Forcing approach is preferred
when the potential inter-user interference is higher than the noise power.
There is a huge performance gain by eliminating the interference.

In the previous method, beams are sent not directly to the users but with some
deviation, which actually reduces the payload. Thirdly, the following
modification, Regularized Zero-Forcing, corrects the beams, which allows some
inter-user interference and significantly increases the payload [21]:

$\bm W_{RZF}=\mu\widetilde{\bm V}^{\mathrm{H}}(\widetilde{\bm V}\widetilde{\bm
V}^{\mathrm{H}}+\bm R)^{-1}\bm P\in\mathbb{C}^{T\times L}.$ (11)

Where diagonal matrix $\bm R\in\mathbb{R}^{L\times L}$ denotes symbol
regularization and matrix $\bm P\in\mathbb{R}^{L\times L}$ and constant
$\mu>0$ as the previous ones. As the baseline, we use a special form [15] of
the regularization matrix $\bm R=\frac{\sigma^{2}L}{P}\bm I$. It is the most
common precoding in real practice.

Algorithm $\bm W_{RZF}$ is a precoding with _scalar regularisation_. Finally,
Adaptive Regularized Zero-Forcing takes into account effective noise of the
system. It was proposed Adaptive RZF (ARZF) algorithm [3] with _diagonal
matrix regularization_ :

$\bm W_{ARZF}=\mu\widetilde{\bm V}^{\mathrm{H}}(\widetilde{\bm
V}\widetilde{\bm V}^{\mathrm{H}}+\lambda\bm S^{-2})^{-1}\bm
P,\quad\lambda=\frac{\sigma^{2}L}{P}.$ (12)

###  3.2 Proposed Optimization Solution

Let us introduce the following parametrization. This method explicitly
constrains antenna rows using (sub-)differentiable projection on a ball as a
part of the Conjugate Spectral Efficiency function: —l— W ∈C ^T ×LS(proj_P,
T(W)), where $\mathcal{S}(\bm W)$ is defined assuming CD $\mathcal{S}(\bm
W)=\textrm{SE}^{C}(\bm W,\widetilde{\bm V},\bm S,\sigma,P)$ or assuming MMSE-
IRC detection $\mathcal{S}(\bm W)=\textrm{SE}(\bm W,\bm H,\bm G^{IRC}(\bm
W),\sigma,P)$ and

$\text{proj}_{P,T}(\bm W)=\begin{cases}\bm w^{m},&\|\bm
w^{m}\|^{2}\leq\frac{P}{T}\\\ \frac{\bm w^{m}}{\|\bm
w^{m}\|}\sqrt{\frac{P}{T}},&\|\bm w^{m}\|^{2}>\frac{P}{T},\quad\forall
m=1\dotsc T\end{cases}$

with $\bm w^{m}\in\mathbb{C}^{L}$ in the $m$-th row of the precoding matrix.

######  Remark 1.

We also add maximization of the target SE function for reference.

######  Remark 2.

As the starting point, we are using Regularized $\bm W_{RZF}$ or Adaptive $\bm
W_{ARZF}$.

The key feature of the proposed method is that the optimization process, on
the one hand, always remains within the feasible range and at the same time
can be solved using unconstrained optimization methods, such as quasi-Newton
methods. In contrast to the projection gradient method, the iteration solution
always remains in a feasible area and does not need an explicit projection on
the boundary. Moreover, in the particular case the projection-based method
doesn’t converge appropriately. And the same for interior point method: its
basic realisation doesn’t give an appropriate solution for our task. An
advantage of the proposed method is that it seeks the solution on the boundary
(and we know that the optimum is on the boundary). This search within the such
complicated edge is possible using the Quasi-Newton L-BFGS method [26].

Input: Channel singular vectors $\widetilde{\bm V}$,

channel singular values $\bm S$, station power $P$, noise $\sigma^{2}$,
iterations $N$,

Tolerance grad $\varepsilon_{g}$, termination tolerance on first order
optimality (default: 1e-5),

Tolerance change $\varepsilon_{c}$, termination tolerance on function value
and parameter changes (default: 1e-9).

Initialize precoding matrix $\bm W\leftarrow\bm W_{0}$;

for  _$t=1$ to $N$_ do

if  _True conditions on $\varepsilon_{g}$ or $\varepsilon_{c}$_ then

return _$\text{proj}_{P,T}(\bm W)$_

end if

Calculate the gradient: $\nabla_{\bm W}\mathcal{S}(\text{proj}_{P,T}(\bm W))$;

Find the optimal direction recursively: $\bm D=\bm D(\nabla_{\bm
W}\mathcal{S}(\text{proj}_{P,T}(\bm W)))$;

Find the optimal step length
$\alpha=\arg\max\limits_{\alpha}\mathcal{S}(\text{proj}_{P,T}(\bm W+\alpha\bm
D))$;

Make the optimization step: $\bm W\leftarrow\bm W+\alpha\bm D$;

end for

return _$\text{proj}_{P,T}(\bm W)$_

Algorithm 1 On the optimal precoding matrix and Quasi-Newton L-BFGS [26]

We consider $\mathcal{S}(\bm W)$ assuming CD or assuming MMSE-IRC detection,
and also initial precoding matrix in the form of RZF or ARZF. By this way, we
compare four options:

  1. 1.

QN CD RZF

  2. 2.

QN CD ARZF

  3. 3.

QN IRC RZF

  4. 4.

QN IRC ARZF

###  3.3 Relative Computational Complexity

The complexity of our algorithm is no more than $\sim 3Nx$ times higher than
the baseline RZF. On each iteration of QNS we have to compute gradient
$\nabla\textrm{SE}^{C}(\bm W)$ of $x$ complexity and adjust step length
$\alpha$, which takes roughly two calls of the function $\textrm{SE}^{C}(\bm
W)$ of $2x$ total complexity (see the Table 1).

RZF | $\textrm{SE}^{C}(\bm W)$ | $\nabla\textrm{SE}^{C}(\bm W)$ | N-iterations of QNS
---|---|---|---
$x$ | $x$ | $x$ | $\sim 3Nx$
Table 1: Relative computational complexity.

###  3.4 Alternative Optimization Solution

For comparison, let us consider the classic softmax parametrization of the
precoding matrix (this is an option of interior point approach). It turned out
that such parametrization leads to the same solution as the differentiable
projection method, but the latter converges at least twice as fast as softmax.
We believe that this is due to the fact that most constraints should be up to
the boundary, and softmax parameterization arguments in this sense should have
the value of infinity. So, we marked the softmax method as an alternative
method and brought it here for a better understanding of the solvable problem.
—l— W ∈C ^T ×LSE^C(W, ~V, S, σ, P)

$\bm W=\\{w_{ij}\\}_{i,j=1}^{T,L},\quad w_{ij}=\rho_{ij}\exp(\bm
i\phi_{ij}),\quad\rho_{ij}\in\mathbb{R}_{+},\quad\phi_{ij}\in[0,2\pi)$
$\rho_{ij}^{2}=\frac{\exp(\theta_{ij})}{\sum_{k=1}^{L}\exp(\theta_{ik})}\frac{\sigma(\alpha_{i})P}{T},\quad\phi_{ij}=2\pi\sigma(\eta_{ij})$

Where $\sigma(\cdot)$ means a sigmoid function, and
$\\{\theta_{ij},\eta_{ij},\alpha_{i}\\}$ are free parameters that can take any
real values. Thus, in another way, we can reduce the problem of optimizing
precoding with constraints to the problem without constraints.

##  4 Numerical Experiments

The main scenario is Urban Non-Line-of-Sight 3GPP_38.901_RMa_NLOS [6] (Fig.
5). We model users in the urban landscape. We are testing the proposed
approach in several scenarios generated using Quadriga [9] — open-source
software for generating realistic radio channel impulse response. For each
scenario, we generate 40 different channel matrices
$H\in\mathbb{C}^{n_{users}\times 4\times 64}$. The carrier frequency for each
channel matrix is selected randomly over the bandwidth. User selection is
described in the next section. The base station antenna array forms a grid
with $8$ placeholders along the $y$ axis and $4$ placeholders along the $y$
axis. The receiver antenna array consists of two placeholders along the $x$
axis. Each placeholder contains two cross-polarized antennas. An interested
reader can find detailed hyperparameters for antenna models and generation
processes in table 3. All unlisted Quadriga parameters are those set by
default.We describe the generation process in detail in our work [3].

For each scenario we generate 40 different channels: $\bm
H\in\mathbb{C}^{K\times R_{k}\times T}$:

  * •

$T=64$ base station antennas,

  * •

$K=8$ users,

  * •

$R_{k}=4$ user antennas,

  * •

$L_{k}=2$ user layers.

Figure 3: Urban NLOS 8 Users. The graph shows IRC Spectral Efficiency of the
different precodings.

###  4.1 Qualitative Results

For each generated channel seed, we fix station power $P=1$ and select system
noise $\sigma^{2}$ for the specified SUSINR (4) average value. We report the
hyper-parameters for Quadriga in the Table 3.

The results of the comparison are in Figure 3 and Table 2. In all cases, the
final SE function is calculated using MMSE-IRC detection (5), geometric mean
effective SINR (2) and Shannon’s capacity (3).

To begin with, the experimental results in Figure 3 and Table 2 show that all
gradient-based methods starting with ”QN,” which stands for ”Quasi-Newton”
(Algorithm (1)), outperform heuristic algorithms such as MRT (9), ZF (10), RZF
(11), and ARZF (12).

The second keyword ”CD” and ”IRC” stands for the assumed detection algorithm:
Conjugate Detection (7) and MMSE-IRC (3) functions, respectively, and their
gradients. It is important to observe that the performance of the obtained
method is influenced by the target optimization function. Remember that the
MMSE-IRC function was used to calculate the final quality measure. As a
result, the ”IRC” method, which optimises the target function directly, has a
higher final quality measure.

On the other hand, the ”CD” method, on the other hand, converges faster
through iterations and each iteration is computationally simpler than IRC. On
the Fig. 4 one can find the quality of the proposed gradient methods by
iterations. We noticed that modification of algorithm that uses approximated
”CD” optimization function gives better results than ”IRC” target function in
the first iterations. At higher iterations the ”IRC” gradient method superiors
the ”CD” method.

The third keyword in the title of gradient methods, ”RZF” (11) and ”ARZF”
(12), refers to the method’s starting point. The method will be more
convergent if the starting point is better. Methods based on the ”ARZF”
starting point outperform methods based on the ”RZF” point in this way.
Remember that our issue isn’t concave. As a result, we only find the local
maximum, and a good starting point is critical for our method to work.

Figure 4: The graph shows IRC Spectral Efficiency on the algorithm iterations.

##  5 Conclusion

The paper studies the multi-user precoding problem as a non-convex
optimization problem for wireless MIMO systems. In this study we approximate
the Spectral Efficiency objective function using a differentiable projection-
based method, where gradients are taken both for the function and the
projection, which tends to very fast convergence of the proposed method. We
then reduce the precoding problem to an unconstrained optimization task and
solve it by the Quasi-Newton L-BFGS iterative procedure. Finally, all our
research is based on realistic power limits per antenna and multiple antenna
users, which in itself is a standalone problem for finding the right approach
to solving. All investigated algorithms were studied in massive experiments
using Quadriga software. The proposed method shows monotonic improvement over
heuristic methods with reasonable computation time. It has a simple
implementation and can be a good reference for other heuristic algorithms in
the MIMO field.

## Acknowledges

Authors are very grateful for D. Minenkov for fruitful discussions.

## Funding

The work is funded by Huawei Technologies.

## References

  * [1] B. Bandemer, M. Haardt, and S. Visuri, _Linear MMSE multi-user MIMO downlink precoding for users with multiple antennas_ , in _2006 IEEE 17th International Symposium on Personal, Indoor and Mobile Radio Communications_. IEEE, 2006, pp. 1–5.
  * [2] E. Björnson, J. Hoydis, and L. Sanguinetti, _Massive mimo networks: Spectral, energy, and hardware efficiency_ , Foundations and Trends in Signal Processing 11 (2017), pp. 154–655.
  * [3] E. Bobrov, B. Chinyaev, V. Kuznetsov, H. Lu, D. Minenkov, S. Troshin, D. Yudakov, and D. Zaev, _Adaptive regularized zero-forcing beamforming in massive mimo with multi-antenna users_ , arXiv preprint arXiv:2107.00853 (2021).
  * [4] E. Bobrov, B. Chinyaev, V. Kuznetsov, D. Minenkov, and D. Yudakov, _Power allocation algorithms for massive mimo systems with multi-antenna users_ (2022).
  * [5] E. Bobrov, D. Kropotov, and H. Lu, _Massive mimo adaptive modulation and coding using online deep learning algorithm_ , IEEE Communications Letters (2021).
  * [6] F. Bohagen, P. Orten, and G. Oien, _Construction and capacity analysis of high-rank line-of-sight MIMO channels_ , in _IEEE Wireless Communications and Networking Conference, 2005_ , Vol. 1. IEEE, 2005, pp. 432–437.
  * [7] J.C. Chen, _Gradient projection-based alternating minimization algorithm for designing hybrid beamforming in millimeter-wave mimo systems_ , IEEE Communications Letters 23 (2018), pp. 112–115.
  * [8] H.H. Dam and A. Cantoni, _Interior point method for optimum zero-forcing beamforming with per-antenna power constraints and optimal step size_ , Signal processing 106 (2015), pp. 10–14.
  * [9] S. Jaeckel, L. Raschkowski, K. Börner, and L. Thiele, _Quadriga: A 3-d multi-cell channel model with time evolution for enabling virtual field trials_ , IEEE Transactions on Antennas and Propagation 62 (2014), pp. 3242–3256.
  * [10] M. Joham, W. Utschick, and J.A. Nossek, _Linear transmit processing in mimo communications systems_ , IEEE Transactions on signal Processing 53 (2005), pp. 2700–2712.
  * [11] M. Joham, W. Utschick, and J.A. Nossek, _Linear transmit processing in mimo communications systems_ , IEEE Transactions on signal Processing 53 (2005), pp. 2700–2712.
  * [12] T.K. Lo, _Maximum ratio transmission_ , in _1999 IEEE international conference on communications (Cat. No. 99CH36311)_ , Vol. 2. IEEE, 1999, pp. 1310–1314.
  * [13] A.H. Mehana and A. Nosratinia, _Diversity of mmse mimo receivers_ , IEEE Transactions on information theory 58 (2012), pp. 6788–6805.
  * [14] L.D. Nguyen, H.D. Tuan, T.Q. Duong, and H.V. Poor, _Multi-user regularized zero-forcing beamforming_ , IEEE Transactions on Signal Processing 67 (2019), pp. 2839–2853.
  * [15] C.B. Peel, B.M. Hochwald, and A.L. Swindlehurst, _A vector-perturbation technique for near-capacity multiantenna multiuser communication-part i: channel inversion and regularization_ , IEEE Transactions on Communications 53 (2005), pp. 195–202.
  * [16] B. Ren, Y. Wang, S. Sun, Y. Zhang, X. Dai, and K. Niu, _Low-complexity mmse-irc algorithm for uplink massive mimo systems_ , Electronics Letters 53 (2017), pp. 972–974.
  * [17] S. Shi, M. Schubert, and H. Boche, _Downlink mmse transceiver optimization for multiuser mimo systems: Duality and sum-mse minimization_ , IEEE Transactions on Signal Processing 55 (2007), pp. 5436–5446.
  * [18] L. Sun and M.R. McKay, _Eigen-based transceivers for the mimo broadcast channel with semi-orthogonal user selection_ , IEEE Transactions on Signal Processing 58 (2010), pp. 5246–5261.
  * [19] S. Verdú, _Spectral efficiency in the wideband regime_ , IEEE Transactions on Information Theory 48 (2002), pp. 1319–1343.
  * [20] B. Wang, Y. Chang, and D. Yang, _On the sinr in massive mimo networks with mmse receivers_ , IEEE Communications Letters 18 (2014), pp. 1979–1982.
  * [21] Z. Wang and W. Chen, _Regularized zero-forcing for multiantenna broadcast channels with user selection_ , IEEE Wireless Communications Letters 1 (2012), pp. 129–132.
  * [22] A. Wiesel, Y.C. Eldar, and S. Shamai, _Zero-forcing precoding and generalized inverses_ , IEEE Transactions on Signal Processing 56 (2008), pp. 4409–4418.
  * [23] D. Wubben, R. Bohnke, V. Kuhn, and K.D. Kammeyer, _Near-maximum-likelihood detection of MIMO systems using MMSE-based lattice-reduction_ , in _2004 IEEE International Conference on Communications (IEEE Cat. No. 04CH37577)_ , Vol. 2. IEEE, 2004, pp. 798–802.
  * [24] T. Yoo and A. Goldsmith, _On the optimality of multiantenna broadcast scheduling using zero-forcing beamforming_ , IEEE Journal on selected areas in communications 24 (2006), pp. 528–541.
  * [25] W. Yu and T. Lan, _Transmitter optimization for the multi-antenna downlink with per-antenna power constraints_ , IEEE Transactions on signal processing 55 (2007), pp. 2646–2660.
  * [26] C. Zhu, R.H. Byrd, P. Lu, and J. Nocedal, _Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization_ , ACM Transactions on mathematical software (TOMS) 23 (1997), pp. 550–560.

Figure 5: Example of Users for Urban Two Building Setup [3]. Precoding | MRT | ZF | RZF | ARZF | QN | QN | QN | QN
---|---|---|---|---|---|---|---|---
SUSINR |  |  |  |  | CD RZF | CD ARZF | IRC RZF | IRC ARZF
-4 | 0.60 | 0.26 | 0.59 | 1.58 | 1.92 | 1.92 | 2.13 | 2.10
0 | 0.79 | 0.52 | 0.85 | 2.25 | 2.70 | 2.69 | 2.95 | 2.95
4 | 0.98 | 0.95 | 1.24 | 3.02 | 3.59 | 3.64 | 3.91 | 3.94
8 | 1.16 | 1.56 | 1.78 | 3.89 | 4.59 | 4.70 | 4.99 | 5.00
12 | 1.33 | 2.38 | 2.52 | 4.86 | 5.73 | 5.85 | 6.17 | 6.20
16 | 1.51 | 3.39 | 3.46 | 5.90 | 6.92 | 7.11 | 7.42 | 7.47
20 | 1.70 | 4.57 | 4.61 | 7.04 | 8.17 | 8.42 | 8.71 | 8.81
24 | 1.91 | 5.93 | 5.95 | 8.29 | 9.43 | 9.83 | 10.01 | 10.25
28 | 2.11 | 7.41 | 7.42 | 9.59 | 10.67 | 11.28 | 11.25 | 11.68
32 | 2.29 | 9.15 | 9.15 | 11.06 | 11.99 | 12.80 | 12.57 | 13.18
36 | 2.45 | 11.06 | 11.06 | 12.63 | 13.53 | 14.35 | 13.97 | 14.64
40 | 2.56 | 13.14 | 13.14 | 14.31 | 15.28 | 15.90 | 15.45 | 16.08
Table 2: Urban NLOS 8 Users. The table shows IRC Spectral Efficiency of the different precoding algorithms. Parameter | Value
---|---
Base station parameters |
number of base stations | $1$
position, m: (x, y, z) axes | $(0,0,25)$
number of antenna placeholders (y axis) | $8$
number of antenna placeholders (z axis) | $4$
distance between placeholders ($y$ axis) |  $0.5$ wavelength
distance between placeholders ($z$ axis) |  $1.7$ wavelength
antenna model | 3gpp-macro
half-Power in azimuth direction, deg | 60
half-Power in elevation direction, deg | 10
front-to back ratio, dB | 20
total number of antennas | 64
Receiver parameters |
number of placeholders at the receiver ($x$ axis) | $2$
distance between placeholders ($x$ axis) |  $0.5$ wavelength
antenna model | half-wave-dipole
total number of antennas | 4
Quadriga simulation parameters |
central band frequency |  $3.5$ GHz
1 sample per meter (default value) | 1
include delay of the LOS path | 1
disable spherical waves (use_3GPP_baseline) | 1
Quadriga channel builders parameters |
shadow fading sigma | 0
cluster splitting | False
bandwidth |  $100$ MHz
number of subcarriers | 42
Table 3: Quadriga generation parameters.