# Novel Approaches to Spectral Properties of Correlated Electron Materials:
From Generalized Kohn-Sham Theory to Screened Exchange Dynamical Mean Field
Theory

Pascal Delange1<EMAIL_ADDRESS>Steffen Backes1
<EMAIL_ADDRESS>Ambroise van Roekeghem2
<EMAIL_ADDRESS>Leonid Pourovskii1,3
<EMAIL_ADDRESS>Hong Jiang4<EMAIL_ADDRESS>and Silke
Biermann1,3,5<EMAIL_ADDRESS>de Physique
Théorique1Centre de Physique Théorique École Polytechnique École Polytechnique
91128 Palaiseau 91128 Palaiseau France
2CEA France
2CEA LITEN LITEN 17 Rue des Martyrs 17 Rue des Martyrs 38054 Grenoble 38054
Grenoble France
3Collège de France France
3Collège de France 11 place Marcelin Berthelot 11 place Marcelin Berthelot
75005 Paris 75005 Paris France
4College of Chemistry and Molecular Engineering France
4College of Chemistry and Molecular Engineering Peking University Peking
University 100871 Beijing 100871 Beijing China
5European Theoretical Spectroscopy Facility (ETSF) China
5European Theoretical Spectroscopy Facility (ETSF) Europe Europe

###### Abstract

The most intriguing properties of emergent materials are typically
consequences of highly correlated quantum states of their electronic degrees
of freedom. Describing those materials from first principles remains a
challenge for modern condensed matter theory. Here, we review, apply and
discuss novel approaches to spectral properties of correlated electron
materials, assessing current day predictive capabilities of electronic
structure calculations. In particular, we focus on the recent Screened
Exchange Dynamical Mean-Field Theory scheme and its relation to generalized
Kohn-Sham Theory. These concepts are illustrated on the transition metal
pnictide BaCo2As2 and elemental zinc and cadmium.

##  1 Introduction

Technological progress has been intimately related with progress in materials
science since its very early days. The last century has seen the development
of refined capabilities for materials elaboration, characterisation and
control of properties, culminating among others in the unprecedented
possibilities of the digital age: materials properties have become an object
of theoretical simulations, stimulating systematic searches for systems with
desired characteristics. This branch of condensed matter physics is
continuously evolving into a true new pillar of materials science, where the
theoretical assessment of solid state systems becomes a fundamental tool for
materials screening.

Obviously, the success of this program hinges on the degree of predictive
power of modern simulation techniques, which have to allow for performing
calculations without introducing adjustable parameters. This requirement is
even more severe since the most popular branch of ab initio calculations, the
Density Functional Theory (DFT) approach, is restricted to ground state
properties of the solid, while most properties of potential technological
interest stem from excited states: the calculation of any type of transport
phenomenon, for example electric, thermal, thermoelectric or magnetoelectric
transport, of optical or spectroscopic properties, or of magnetic, charge- or
orbital susceptibilites requires the accurate assessment of highly non-trivial
response functions within a finite-temperature description.

It is not a coincidence that the materials with the most exotic features are
also most challenging to get to grips with from a theoretical point of view:
both the complexity of their properties and the difficulty of their
description are in fact a consequence of the very nature of the underlying
electronic structure: typically, one is dealing with transition metal,
lanthanide or actinide compounds, where the electrons in partially filled
$d$\- or $f$-shells strongly interact with each other through electron-
electron Coulomb interactions, leading to highly entangled many-body quantum
states.

Even when using a simplified description of the solid in terms of an effective
lattice model, assessing these quantum correlated states is a tremendous
challenge. If local quantum fluctuations on a given atomic site are dominant,
dynamical mean-field theory (DMFT) [1, 2] yields an accurate description of
the system. In particular, DMFT captures both the strong coupling Mott-
insulating limit and the weakly interacting band limit of the Hubbard model,
and is able to describe the salient spectral features of a correlated metal
with coexisting quasi-particle and Hubbard peaks even in the intermediate
correlation regime. In the multi-orbital case, additional degrees of freedom
can lead to even richer physics with unconventional (e.g. orbital-selective)
behavior [3], or complex ordering phenomena [4, 5, 6, 7, 8, 9].

In order to recover the material-specific character of the calculations, DMFT
has been combined with DFT [10, 11] into the so-called “DFT+DMFT” scheme,
which is nowadays one of the most popular workhorses of electronic structure
theory for correlated electron materials. Successful applications include
transition metals [12, 13, 14], transition metal oxides [4, 15, 16, 17, 18,
19, 20, 21, 22, 23, 24, 25, 26], lanthanide [27, 28] or actinide [29, 30, 31,
32] systems. Comparisons between calculated and measured spectral functions
have sometimes led to impressive agreement[33, 34]. However, despite these
successes the construction of the Hamiltonian used in these calculations is
rather ad hoc and remains a limitation to the predictive power of the
approach. Therefore, the elaboration of more systematic interfaces between
electronic structure and many-body theory has become an active area of
research [35].

In this work, we review a recent scheme combining screened exchange and DMFT
[36, 37, 38]. We discuss its relation to many-body perturbation theory and
generalized Kohn-Sham (KS) Theory and analyse the effects included at the
different levels of the theory. As an illustration we describe its application
to the low-energy spectral properties of the cobalt pnictide BaCo2As2.

How does a theory thas has been designed for strongly correlated materials
reduce to an a priori simpler version in the case of a weakly correlated
system ? For methods based on DFT that include only the static on-site Coulomb
interaction between localized states, like DFT+U[39] or DFT+DMFT, the spectrum
in the limit of vanishing Hubbard interactions reduces trivially to the Kohn-
Sham spectrum of DFT. For more complex interfaces of electronic structure and
many-body theory which we will discuss here, however, it is a non-trivial
question and allows for interesting possibilities to check their consistency.
Here, we will discuss this issue on the example of the transition metals zinc
and cadmium. This also allows us to comment on the challenge of including
states in a wider energy range, identifying a challenging problem on these
seemingly “simple” systems.

The paper is organised as follows: in section 2 we review how Screened
Exchange Dynamical Mean Field Theory derives from the combined many-body
perturbation theory \+ dynamical mean field scheme “GW+DMFT”. In section 3 we
analyse the relation of screened exchange schemes to generalized KS theory.
Section 4 provides an example of the application of such methods to BaCo2As2,
while section 5 discusses the electronic structure of elemental zinc and
cadmium within a simplified approach. Finally, we conclude in section 6.

##  2 From GW+DMFT to Screened Exchange Dynamical Mean Field Theory

Improving the predictive power of methods that combine electronic structure
and many-body theory poses the challenge of properly connecting the two
worlds, without double counting of interactions or screening. At the heart of
this challenge lies the mismatch between the density-based description of DFT
and the Green’s function formalism used at the many-body level, as well as the
difficulty of incorporating the feedback of high-energy screening processes
governed by the unscreened Coulomb interaction onto the low-energy electronic
structure. Conceptually speaking, these difficulties can be avoided by working
on a large energy scale in the continuum with the full long-range Coulomb
interactions, and a Green’s function-based formalism even at the level of the
weakly correlated states. These features are realised within the combined
many-body perturbation theory and dynamical mean field theory scheme “GW+DMFT”
[40, 41, 42, 43, 44, 45, 46, 47, 48, 49]: screening is assessed by the random
phase approximation in the continuum, augmented by a local vertex correction,
while the starting electronic structure for the DMFT calculation can be
roughly interpreted as a “non-local GW” calculation [50].

Screened Exchange Dynamical Mean Field Theory [36, 37, 38] can be understood
as an approximation to this full GW+DMFT scheme. It is based on the recent
observation [36, 50] that within the GW approximation the correction to LDA
can be split into two contributions: a local dynamical self-energy
$\Sigma_{loc}(\omega)$ and a k-dependent but static self-energy
$\Sigma_{nloc}(k)$, which does not contain any local component. If such a
separation was strictly valid in the full energy range, the non-local static
part of the full self-energy would be given by the non-local Hartree-Fock
contribution, since the dynamical part vanishes at high frequency. In many
realistic systems such a decomposition holds to a good approximation in the
low-energy regime that we are interested in, where the static part
$\Sigma_{nloc}(k)$ is quite different from the Fock exchange term. It is
approximately given by a screened exchange self-energy, leading to a
decomposition of the GW self-energy into
$\Sigma_{GW}=[GW(\nu=0)]_{nonloc}+[GW]_{loc}$. Here, the first term is a
screened exchange contribution arising from the screened interaction $W(\nu)$
evaluated at zero frequency $W(\nu=0)$. The second term is the local
projection of the $GW$ self-energy. It is simply given by the GW self-energy
evaluated using a local propagator $G_{loc}$ and the local screened Coulomb
interaction $W_{loc}$. Exactly as in GW+DMFT, Screened Exchange Dynamical Mean
Field Theory replaces this term by a non-perturbative one: it is calculated
from an effective local impurity problem with dynamical interactions. In
current practical applications of Screened Exchange + DMFT the RPA-screened
Coulomb potential $W$ has been replaced by its long-wavelength limit, which
reduces to a simple Yukawa-form [36].

Quite generally the dynamical character of the interactions results in an
additional renormalization $Z_{B}$ of the hopping amplitudes, which can be
understood as an electronic polaron effect: the coupling of the electrons to
plasmonic screening degrees of freedom leads to an effective mass enhancement
corresponding to the hopping reduction, manifesting itself as a narrowing of
the band. This effect can be estimated from the plasmon density of modes as
given by the imaginary part of the frequency dependent interaction $W$. An
explicit expression for $Z_{B}$ has been derived in Ref.
PhysRevLett.109.126408.

##  3 Relation to Generalized Kohn-Sham Theory

Screened Exchange Dynamical Mean Field Theory can also be viewed as a specific
approximation to a spectral density functional theory based on the Generalized
Kohn-Sham (GKS) scheme of Seidl et al.[52]. In GKS theory, alternative choices
for the reference system that are different than the familiar Kohn-Sham system
of DFT are explored. In particular, a generalized Kohn-Sham scheme where the
reference system is a screened exchange Hamiltonian can be constructed. The
main motivation for the inclusion of screened exchange in the literature has
been to improve upon the band gap problem in semiconductors. Indeed, it can be
shown that the screened exchange contribution, which corresponds to a non-
local potential, effectively reintroduces to some degree the derivative
discontinuity that is missing in the pure DFT description based on local
exchange-correlation potentials[53]. Since the derivative discontinuity
corresponds to the discrepancy between the true gap and the Kohn-Sham gap in
exact DFT, a substantial improvement of the theoretical estimate for band gaps
can be expected on physical grounds and has indeed been found. Here, our goal
is somewhat different: motivated by the analysis of the role of screened
exchange in GW+DMFT described above, we would like to connect the Screened
Exchange DMFT scheme introduced above to generalized KS schemes making direct
use of the non-local screened exchange potential.

With this aim in mind, we briefly review the generalized Kohn-Sham
construction in the case an effective Kohn-Sham system including screened
exchange. Hereby, we follow closely Seidl et al.[52], both in notation and
presentation. First, one defines a functional

$\displaystyle
S\left[\Phi\right]=\left\langle\Phi\left|\hat{T}\right|\Phi\right\rangle+U_{H}\left[\\{\phi_{i}\\}\right]+E^{sx}_{x}\left[\\{\phi_{i}\\}\right]$
(1)

that includes, in addition to the familiar kinetic energy term
$\langle\Phi|\hat{T}|\Phi\rangle$ and the Hartree energy
$U_{H}\left[\\{\phi_{i}\\}\right]$ also the screened Fock term

$\displaystyle E^{sx}_{x}\left[\\{\phi_{i}\\}\right]=$
$\displaystyle-\sum\limits_{i<j}^{N}\int
d\boldsymbol{r}d\boldsymbol{r^{\prime}}$
$\displaystyle\times\cfrac{\phi_{i}^{*}(\boldsymbol{r})\phi_{j}^{*}(\boldsymbol{r^{\prime}})e^{-k_{\mathrm{TF}}\left|\boldsymbol{r}-\boldsymbol{r^{\prime}}\right|}\phi_{j}(\boldsymbol{r})\phi_{i}(\boldsymbol{r^{\prime}})}{\left|\boldsymbol{r}-\boldsymbol{r^{\prime}}\right|}$
(2)

Here, $\Phi$ are Slater determinants of single-particle states $\phi_{i}$.
$k_{\mathrm{TF}}$ is the Thomas Fermi wave vector. In order to derive a
functional of the density Seidl et al. define a functional $F^{s}$ via the
minimisation

$\displaystyle
F^{S}\left[\rho\right]=\min_{\Phi\rightarrow\rho(\boldsymbol{r})}S\left[\Phi\right]=\min_{\left\\{\phi_{i}\right\\}\rightarrow\rho(\boldsymbol{r})}S\left[\left\\{\phi_{i}\right\\}\right]$
(3)

Next we define the energy functional

$\displaystyle
E^{S}\left[\\{\phi_{i}\\};v_{eff}\right]=S\left[\\{\phi_{i}\\}\right]+\int
d\boldsymbol{r}v_{eff}(\boldsymbol{r})\rho(\boldsymbol{r})$ (4)

where now the potential $v_{eff}$ does not only include the external potential
$v$ as in usual DFT, but also a contribution by the exchange-correlation part

$\displaystyle v_{eff}=v+v_{xc}^{sx}\left[\rho\right].$ (5)

The additional contribution, the generalized (local) exchange-correlation
potential

$\displaystyle v_{xc}^{sx}=\frac{\partial
E_{xc}^{sx}\left[\rho\right]}{\partial\rho}.$ (6)

is the functional derivative of the generalized (local) exchange-correlation
functional

$\displaystyle
E_{xc}^{sx}\left[\rho\right]=E_{xc}\left[\rho\right]-E_{x}^{sx}\left[\rho\right]+T\left[\rho\right]-T^{sx}\left[\rho\right]$
(7)

which comprises the difference between the exchange-correlation potential of
standard Kohn-Sham DFT and the non-local exchange energy defined above, as
well as the difference between the kinetic energies of the standard and
generalized Kohn-Sham systems. The functional derivative will eventually have
to be evaluated self-consistently at the converged density.

This construction leads to the generalized Kohn-Sham equations

$\displaystyle-\nabla^{2}$
$\displaystyle\phi_{i}(\boldsymbol{r})+v(\boldsymbol{r})\phi_{i}(\boldsymbol{r})+u\left([\rho];\boldsymbol{r}\right)\phi_{i}(\boldsymbol{r})$
$\displaystyle-\int
d\boldsymbol{r^{\prime}}v_{x}^{sx,\mathrm{NL}}(\boldsymbol{r},\boldsymbol{r^{\prime}})\phi_{i}(\boldsymbol{r^{\prime}})+v_{xc}^{sx}\left([\rho];\boldsymbol{r}\right)\phi_{i}(\boldsymbol{r})=\varepsilon_{i}\phi_{i}$
(8)

with the Hartree potential $u$ and the non-local screened Fock potential

$\displaystyle v_{x}^{sx,\mathrm{NL}}(\boldsymbol{r},\boldsymbol{r^{\prime}})$
$\displaystyle=-\sum\limits_{j=1}^{N}\cfrac{\phi_{i}(\boldsymbol{r})e^{-k_{\mathrm{TF}}\left|\boldsymbol{r}-\boldsymbol{r^{\prime}}\right|}\phi_{j}^{*}(\boldsymbol{r^{\prime}})}{\left|\boldsymbol{r}-\boldsymbol{r^{\prime}}\right|}.$
(9)

and the effective (local) generalized Kohn Sham potential $v_{xc}^{sx}$
defined above. The generalized Kohn-Sham equations have the form

$\displaystyle\hat{\mathcal{O}}\left[\\{\phi_{i}\\}\right]\phi_{j}+\hat{v}_{eff}\phi_{j}=\epsilon_{j}\phi_{j}$
(10)

where $\hat{\mathcal{O}}$ is a non-local operator, generalizing the standard
Kohn-Sham operator consisting solely of kinetic energy and Hartree potential.

The ground state energy for a system in the external potential $v$ is then
given by the expression

$\displaystyle E_{{\rm SEx-
DFT}}\left[v\right]=F^{S}\left[\rho_{0}^{S}\left[v_{eff}\right]\right]+E_{xc}^{sx}\left[\rho_{0}^{S}\left[v_{eff}\right]\right]$
$\displaystyle+\int drv(r)\rho_{0}^{S}\left[v_{eff}\right]$ (11)

Now, the relation to Screened Exchange DMFT is becoming clear: one may
construct a spectral density functional in the same spirit as in DFT+DMFT
[54], but starting from the generalized KS functional. In amadon:066402, the
expression for the total energy within the standard DFT+DMFT case was derived
to be

$E=E_{{\rm DFT}}-\sum_{l}\epsilon^{{\rm KS}}_{l}+\langle H_{{\rm
KS}}\rangle+\langle(H_{int}-H_{dc})\rangle$ (12)

where $\sum_{l}\epsilon^{{\rm KS}}_{l}$ is the sum of the occupied Kohn-Sham
eigenvalues, $\langle H_{{\rm KS}}\rangle=\mbox{tr}\,[H_{KS}\hat{G}]$, and
$H_{int}$ and $H_{dc}$ denote the local interaction part of the Hamiltonian
and the corresponding double counting term, respectively.

Instead of using the usual Kohn-Sham Hamiltonian for the construction of the
one-body part, Screened Exchange DMFT relies on the generalized Kohn-Sham
reference system that includes the screened exchange potential. The
generalization of (12) to the present case thus replaces the Kohn-Sham
Hamiltonian $H_{{\rm KS}}$ in the expression for the energy by its non-local
form, keeping track of the effective potential part:

$\displaystyle E=E_{{\rm SEx-DFT}}-\sum_{l}\epsilon^{{\rm SEx-
KS}}_{l}+\langle\hat{\mathcal{O}}+\hat{v}_{eff}\rangle$
$\displaystyle+\langle(H_{int}(V_{ee},\lambda_{s},\omega_{s})-H_{dc})\rangle$
(13)

Furthermore, the local interaction term is taken in the more general form of a
dynamical interaction, thus corresponding to a local Hubbard term with
unscreened interactions and local Einstein plasmons of energy $\omega_{s}$
coupling to the electrons via coupling strength $\lambda_{s}$.

Figure 1: (Color online) Kohn-Sham band structure of BaCo2As2 within DFT-LDA
(red lines) and the screened exchange approximation (black lines).

This concludes our description of the generalized Kohn-Sham interpretation of
Screened Exchange DMFT, resulting in particular in an energy functional
expression. However, in the practical calculations presented in the following
we do not minimise the full energy expression as given above, but rather work
at the converged DFT density and then investigate spectral properties using
the Screened Exchange DMFT formalism. This amounts to a one-shot Screened
Exchange-DFT+DMFT calculation that uses the DFT density as a starting point.
The advantage of such an approach is obvious: Numerically, this procedure
allows us to avoid the expensive evaluation of non-local exchange terms within
the self-consistency cycle of GKS theory. Moreover, as is well-known, while
severe deviations of the true spectrum from the Kohn-Sham spectrum are quite
common, the ground state density obtained even from approximate DFT
functionals is often a good representation of the true one. In the case of the
exact DFT functional, our approach would also lead to the exact ground state
density and energy, with additional improvements of the spectrum over standard
Kohn-Sham DFT.

Figure 2: (Color online) Kohn-Sham band structure of BaCo2As2 within DFT-LDA,
projected on the $d_{x^{2}-y^{2}}$ orbital.

##  4 Results on BaCo2As2

As a first illustration of Screened Exchange DMFT, we review calculations on
BaCo2As2, which is the fully Co-substituted representative of the so-called
“122” family of the iron-based superconductors, isostructural to the
prototypical parent compound BaFe2As2. The Fe $\rightarrow$ Co substitution
has however important consequences: the nominal filling of the $3d$ states
changes from d6 to d7, which strongly reduces the degree of Coulomb
correlations, making BaCo2As2 a moderately correlated compound[56, 36].
Indeed, the power-law deviations from Fermi liquid behavior above an extremely
low coherence temperature discussed in the Fe-based compounds[57] are a
consequence of the d6 configuration and strong intra-atomic exchange
interactions. This is no longer the case in the cobalt pnictides, where angle-
resolved photoemission spectroscopy (ARPES) identifies well-defined and long-
lived quasiparticle excitations with relatively weak mass renormalization.
Nevertheless, the DFT-LDA derived Fermi surface differs from experiment[56,
58, 36, 59]. Therefore, BaCo2As2 provides an ideal testing ground for new
approaches to spectroscopic properties. We do not discuss here the details of
yet another interesting question, which is the absence of ferromagnetism
despite a high value of the DFT density of states at the Fermi level,
suggestive of Stoner ferromagnetism, but refer the reader to Ref.
PhysRevLett.113.266403, where the solution to this puzzle was discussed in
detail.

Figure 3: (Color online) Angle-resolved photoemission spectrum of BaCo2As2.
Adapted from Ref. PhysRevLett.113.266403.  Figure 4: (Color online) k-resolved
spectral function of BaCo2As2 within LDA+DMFT. Adapted from Ref.
PhysRevLett.113.266403. Figure 5: (Color online) k-resolved spectral function
of BaCo2As2 within Screened Exchange DMFT (with dynamical interactions).
Adapted from Ref. PhysRevLett.113.266403. Figure 6: (Color online) k-resolved
spectral function of BaCo2As2 within Screened Exchange DMFT (with dynamical
interactions), projected on the $d_{x^{2}-y^{2}}$ orbital. Adapted from Ref.
PhysRevLett.113.266403.

Fig. 1 shows the DFT band structure of BaCo2As2, in comparison to a screened
exchange calculation. As in the iron-based pnictides, the dominantly
$3d$-derived states are located around the Fermi level; in this case in a
window of about -3 eV to 2 eV. As compared to the parent iron pnictides with
$3d^{6}$ configuration of the Fe shell, the Fermi surface topology is modified
due to the larger $3d^{7}$ filling. The hole pocket at the $\Gamma$ point that
is present in most Fe-based pnictide compounds is pushed below the Fermi
level, as well as the band forming the electron pocket at $M$, which is now
fully filled. In standard DFT-LDA (see Fig. 2), a characteristic flat band of
dominant $x^{2}-y^{2}$ character lies directly on the Fermi level around the M
point, giving rise to a huge peak in the density of states. The close
proximity of this band to the Fermi level will render its energy highly
sensitive to the details of the calculation and thus provides a perfect
benchmark for improved computational techniques.

Comparison to angle-resolved photoemission experiments[36], as reproduced in
Fig. 3, reveals that the overall bandwidth of the LDA band structure is too
wide by roughly a factor of 1.5. In a DFT+DMFT calculation, as shown in Fig.
4, this band renormalisation is reproduced, giving an overall occupied
bandwidth of about 1.5 eV. The fine details of the Fermi surface, and in
particular the position of the flat $x^{2}-y^{2}$ band are however not well
described. For a detailed comparison we refer the reader to Ref.
PhysRevLett.113.266403.

Including screened exchange in the form of a Yukawa potential on top of DFT,
as shown in Fig. 1 widens the band by a considerable amount as compared to
DFT-LDA. Even more striking are the modifications at the Fermi level: the
$x^{2}-y^{2}$ band has been shifted above the Fermi level, with an energy at
the $M$ point of about 0.15 eV above $E_{F}$. While the lowering of the
filling of this band improves the agreement between theory and experiment, the
shift is too large to reproduce the experimental Fermi surface. However, when
applying DMFT with dynamical interactions on top of the screened exchange
Hamiltonian, equivalent to the one-shot Screened Exchange DMFT procedure
described above, this band is renormalized by the electronic interactions and
ends up again close (but above!) the Fermi level. Its energy is slightly
higher than in DFT-LDA and is now in excellent agreement with the experimental
Fermi surface. A more detailed analysis presented in Ref. [36] reveals that
also the higher energy features such as the bands within the range of up to 2
eV are well reproduced in this one-shot Screened Exchange DMFT approach.

The comparison of the band structures and spectral functions within the
different computational schemes illustrates the effects of the different terms
in an instructive way: Improving on the description of exchange by replacing
the local exchange as contained in DFT-LDA by a non-local screened Fock
exchange widens the bandwidth and strongly “overcorrects” the Fermi surface.
Improving at the same time on the correlation part by applying DMFT with
frequency-dependent interactions then leads to a partial cancellation of the
band widening, giving an overall bandwidth surprisingly close to the LDA+DMFT
one. This suggests that the good description of the overall bandwidth in
DFT+DMFT is the result of an error cancellation between the local
approximation to exchange and partial neglect of correlations. However, the
Fermi surface is strongly modified by the non-local corrections, leading to a
significant improvement over LDA+DMFT and resulting in good agreement with
experiment.

These findings might suggest that possible inconsistencies between theoretical
and experimental Fermi surfaces that are observed in many pnictide materials
could be corrected by adding simple non-local self-energy corrections stemming
from screened exchange effects. Whether this treatment can fully account for
the “red-blue shift” of Fermi surface pockets found in the experimental-
theoretical comparisons[60] is a most interesting topic for future research.
The current example may give rise to optimism.

##  5 How does Screened Exchange Dynamical Mean Field Theory behave for weakly
correlated materials?

Figure 7: (Color online) The band structure of elemental Zn calculated within
DFT. The orbital character is indicated by the intensity of the different
colors.  Figure 8: (Color online) The band structure of elemental Cd
calculated within DFT. The orbital character is indicated by the intensity of
the different colors.

Hamiltonians built as combinations of a DFT part and local Hubbard-type
interaction terms trivially reduce to the DFT Kohn-Sham electronic structure
when assuming that in weakly correlated materials the effective local
interactions become negligible. The question of the recovery of the weakly
interacting limit is, however, more interesting in the case of Screened
Exchange DMFT. While the static part of the effective local interaction may be
assumed to lose its importance, band widening by the replacement of the DFT
exchange correlation potential by the non-local exchange-correlation GW self-
energy persists. On the other hand, plasmonic effects are also present in
weakly correlated materials and continue to renormalize the low-energy band
structure through electron-plasmon coupling. This raises the question of what
the resulting spectra for weakly correlated materials look like in screened
exchange + DMFT.

In Ref. 0295-5075-108-5-57003, this question has been studied for early
transition metal perovskites, where it was found that the band widening effect
induced by non-local exchange and the electronic polaron effect counteract
each other and tend to approximately cancel, thus resulting in a low-energy
electronic structure close to the DFT Kohn-Sham band structure as long as
static Hubbard interactions are disregarded. Here, we address this question in
the case of the seemingly “simple” transition metals zinc and cadmium.

Both elements nominally display a $d^{10}$ configuration, with fully occupied
$3d$ orbitals in the case of Zn and $4d$ in the case of Cd, the dominantly
$d$-derived bands being located several eV below the Fermi level. In Fig. 7
and Fig. 8 we show the band structure calculated within DFT for both
materials. Here and in the following we use the experimental crystal
structure. DFT puts the occupied $d$ states at around -8 eV in Zn and -9 eV in
Cd. The conducting states of this transition metal are formed by dispersive
4$s$(5$s$) in Zn(Cd) states around the Fermi level, that hybridize with the
$p$-manifold.

These facts raise the immediate expectation of negligibly small correlation
effects on the occupied $d$ shells. An effective Hubbard interaction
calculated for the $d$-manifold within the constrained random phase
approximation coincides with the fully screened interaction, since as a
consequence of the complete filling of the d-shell there are no intra-d
transitions to be cut out as opposed to an open shell system, where
transitions inside the shell contribute to screening effects.

Figure 9: (Color online) The fully screened effective local Hubbard
interaction on the $3d$ manifold for Zn. Figure 10: (Color online) The fully
screened effective local Hubbard interaction on the $4d$ manifold for Cd.

Figs. 9 and 10 display the local component of this fully screened interaction
projected on the $d$-manifold. The low-frequency limit approaches a value of
4.8 eV and 3.1 eV for Zn and Cd respectively. Even though this value is
similar to their oxides, which are open shell systems, where correlations on
the $d$ states are significant, the high binding energy of these states far
away from the Fermi level effectively prevents dynamical fluctuations. This
suggests that screened exchange + DMFT should in fact reduce to screened
exchange renormalized by the bosonic factor $Z_{B}$ discussed above.
Nevertheless, this does not mean that static effects of the interaction are
properly treated in DFT. Even if this were the case, there is no reason that
the DFT Kohn-Sham spectrum, being derived from an effective non-interacting
system, provides an accurate description of the experimental situation.

Figure 11: (Color online) The Density of States of elemental Zn calculated
within Density Functional Theory (DFT, black solid line) and Hartree-Fock (HF,
red solid line) in comparison with Photoemission experiments [61, 62] (dashed
line, symbols). Density Functional Theory (DFT) calculations underestimate the
binding energy of the Zn $3d$ states, while the HF overestimates the binding
energy significantly (see explanation in the text). Figure 12: (Color online)
The Density of States of elemental Cd calculated within Density Functional
Theory (DFT, black solid line) and Hartree-Fock (HF, red solid line) in
comparison with Photoemission experiments [61] (dashed line, symbols). Density
Functional Theory (DFT) calculations underestimate the binding energy of the
Cd $3d$ states, while the HF overestimates the binding energy significantly
(see explanation in the text).

Figs. 11 and 12 compare the experimental photoemission spectra [61, 62] from
the literature to the density of states (DOS) derived from DFT and Hartree-
Fock (HF) theory. The resulting discrepancy in terms of an underestimation of
the binding energy of the $d$ states in DFT of several eV had been noted in
the literature before[63, 64, 65]: Norman et a.l [63] discussed it in terms of
a self-interaction error, proposing a correction in terms of an approximate
substraction of self-interaction contributions contained in DFT [66]. Hartree-
Fock calculations are, on the one hand, self-interaction free, but on the
other hand – due to the absence of screening – widen all bands and place the
$d$-bands far too low in energy, as can be seen in Figs. 11 and 12.

Since local dynamical correlations can assumed to be small as discussed
before, an improved treatment of the screened interaction and the self-
interaction correction at the same time is likely to improve the shortcomings
of Hartree-Fock, resp. DFT. Here we will discuss the two possible extensions
of Screened Exchange plus a bosonic renormalization factor $Z_{B}$ and the GW
approximation. The computationally cheaper option of screened exchange is
including only static exchange contributions with a Yukawa-type interaction
potential, and an effective renormalization $Z_{B}$ which originates from the
spectral weight transfer to plasmonic excitations. GW is computationally more
demanding, but has the advantage of treating the dynamical part of the
screening and correlation. Even though both methods include a self-interaction
correction, the self-interaction contained in the Hartree term is not
completely cancelled since the exchange contributions are derived from a
screened interaction and not the bare one, as opposed to the Hartree term.

In Figs. 13 and 14 we show comparisons of the DOS of zinc and cadmium,
calculated within DFT, screened Exchange(+$Z_{B}$) and GW[67] to photoemission
spectra. In the GW calculation we used $7\times 7\times 3$ k-points and 5
additional high-energy local orbitals. Interestingly, while in both systems
the GW approximation provides a significant correction of the DFT Kohn-Sham
spectrum in the right direction, but still underestimates the binding energy,
the screened exchange scheme places the d-states too low in energy for Zn
while providing a slightly better estimate than GW in Cd. The addition of the
bosonic renormalization factor $Z_{B}$ merely renormalizes the d bandwidth but
keeps the average level position constant.

Figure 13: (Color online) The Density of States of elemental Zn calculated
within different theoretical methods (solid lines), in comparison with
Photoemission experiments [61, 62] (dashed line, symbols). Density Functional
Theory (DFT) calculations significantly underestimate the binding energy of
the Zn $3d$ states, while the GW approximation obtains a much better
agreement. Screened exchange overestimates the binding energy significantly.
Figure 14: (Color online) The Density of States of elemental Cd calculated
within different theoretical methods (solid lines), in comparison with
Photoemission experiments [61] (dashed line). DFT calculations significantly
underestimate the binding energy of the Cd $4d$ states, while the GW
approximation and also screened exchange obtain a much better agreement.

This raises the interesting question of which effects are missing in screened
exchange and GW? The incomplete cancellation of the self-interaction in both
approaches is expected to lead to an overall underestimation of the binding
energy, since the additional unphysical interaction increases the energy of
the d states. A more accurate estimate of this term would lead to an
improvement of GW in both systems, but an even larger error of screened
exchange in Zn. Another effect neglected in screened exchange is the Coulomb
hole contribution: this term, discussed by Hedin as part of the “Coulomb hole
screened exchange (COHSEX)” approximation translates the fact that the
presence of an electron at a position $r$ pushes away charge at $r$ (in the
language of a lattice model, the charge-charge correlation function exhibits a
reduction of the double occupancy), and the interaction of this effective
positive charge with the electron presents an energy gain expressed in the
form of an interaction of the electron with a “Coulomb hole”. This term,
contained in GW but not in screened exchange, also increases the binding
energy.

This leads to the overall picture that screened exchange with the inclusion of
the static corrections just discussed has a tendency to overestimate the
binding energy in general, while GW underestimates it. The obvious difference
between the two methods is the dynamical treatment of the screened
interaction, which is treated more appropriately in GW, but it is not clear a
priori whether the static approximation of the screened interaction or the
approximated form of the screened interaction in terms of a Yukawa potential
gives rise to the difference between screened exchange and GW.

The GW description of zinc and cadmium is close to the experiment. The
remaining discrepancy to experiment is likely explained by remaining self-
interaction contributions and/or missing self-consistency. Self-consistency
(or quasiparticle self-consistency) has been investigated in the homogeneous
electron gas [68] and various solid state systems [69]. These questions are
left for future work.

##  6 Conclusions

In this work we have reviewed and applied existing as well as novel approaches
to obtain spectral properties of correlated electron materials. Guided by the
need of a proper treatment of the long-range Coulomb interaction and non-local
exchange effects we presented a lightweight version of the general GW+DMFT
approach, the so-called Screened Exchange DMFT. It can be derived as a
simplication to GW+DMFT in terms of a generalized screened exchange DFT scheme
where local interactions are treated by dynamical DMFT.

Analysis and the application of a simplifed form of this scheme to BaCo2As2
indeed showed that non-local exchange and electronic screening lead to
significant corrections to the electronic spectrum which are necessary to
obtain a proper description of the experimental observations.

Furthermore, we discussed the case of the elemental transition metals Zn and
Cd, where strong local correlations are unimportant but the position of the
occupied d manifold is very sensitive to a proper treatment of screened
exchange effects.

SB and AvR thank their collaborators of Refs.36 and 37 for the fruitful
collaborations and discussions. This work was supported by IDRIS/GENCI Orsay
under project t2017091393, the European Research Council under its
Consolidator Grant scheme (project 617196), and the ECOS-Sud grant A13E04. HJ
acknowledges financial support from National Natural Science Foundation
(21373017, 21621061) and the National Basic Research Program of China
(2013CB933400).

## References

  * [1] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg: Reviews of Modern Physics 68 (1996) 13.
  * [2] G. Kotliar and D. Vollhardt: Physics Today March 2004 (2004) 53.
  * [3] S. Biermann, L. de’ Medici, and A. Georges: Phys. Rev. Lett. 95 (2005) 206401.
  * [4] E. Pavarini, S. Biermann, A. Poteryaev, A. I. Lichtenstein, A. Georges, and O. K. Andersen: Phys. Rev. Lett. 92 (2004) 176403.
  * [5] C. Martins, M. Aichhorn, L. Vaugier, and S. Biermann: Phys. Rev. Lett. 107 (2011) 266404.
  * [6] R. Arita, J. Kuneš, A. V. Kozhevnikov, A. G. Eguiluz, and M. Imada: Phys. Rev. Lett. 108 (2012) 086403.
  * [7] K.-H. Ahn, K.-W. Lee, and J. Kuneš: Journal of Physics: Condensed Matter 27 (2015) 085602.
  * [8] J. Kuneš and V. Křápek: Phys. Rev. Lett. 106 (2011) 256401.
  * [9] J. Kuneš: Phys. Rev. B 90 (2014) 235140\.
  * [10] V. Anisimov, A. Poteryaev, M. Korotin, A. Anokhin, and G. Kotliar: Journal of Physics: Condensed Matter 9 (1997) 7359.
  * [11] A. Lichtenstein and M. Katsnelson: Physical Review B 57 (1998) 6884\.
  * [12] A. I. Lichtenstein, M. I. Katsnelson, and G. Kotliar: Phys. Rev. Lett. 87 (2001) 067205.
  * [13] P. Delange, T. Ayral, S. I. Simak, M. Ferrero, O. Parcollet, S. Biermann, and L. Pourovskii: Phys. Rev. B 94 (2016) 100102.
  * [14] J. Sánchez-Barriga, J. Minár, J. Braun, A. Varykhalov, V. Boni, I. Di Marco, O. Rader, V. Bellini, F. Manghi, H. Ebert, M. I. Katsnelson, A. I. Lichtenstein, O. Eriksson, W. Eberhardt, H. A. Dürr, and J. Fink: Phys. Rev. B 82 (2010) 104414.
  * [15] S. Biermann, A. Poteryaev, A. I. Lichtenstein, and A. Georges: Phys. Rev. Lett. 94 (2005) 026404.
  * [16] A. I. Poteryaev, J. M. Tomczak, S. Biermann, A. Georges, A. I. Lichtenstein, A. N. Rubtsov, T. Saha-Dasgupta, and O. K. Andersen: Phys. Rev. B 76 (2007) 085127.
  * [17] J. M. Tomczak and S. Biermann: Journal of Physics: Condensed Matter 19 (2007) 365206.
  * [18] J. M. Tomczak and S. Biermann: EPL (Europhysics Letters) 86 (2009) 37004\.
  * [19] J. M. Tomczak and S. Biermann: Journal of Physics: Condensed Matter 21 (2009) 064209.
  * [20] J. M. Tomczak, A. I. Poteryaev, and S. Biermann: Comptes Rendus Physique 10 (2009) 537 .
  * [21] I. A. Nekrasov, N. S. Pavlov, and M. V. Sadovskii: Journal of Experimental and Theoretical Physics 116 (2013) 620.
  * [22] P. Thunström, I. Di Marco, and O. Eriksson: Phys. Rev. Lett. 109 (2012) 186401.
  * [23] I. Leonov, L. Pourovskii, A. Georges, and I. A. Abrikosov: Phys. Rev. B 94 (2016) 155135.
  * [24] V. I. Anisimov, D. E. Kondakov, A. V. Kozhevnikov, I. A. Nekrasov, Z. V. Pchelkina, J. W. Allen, S.-K. Mo, H.-D. Kim, P. Metcalf, S. Suga, A. Sekiyama, G. Keller, I. Leonov, X. Ren, and D. Vollhardt: Phys. Rev. B 71 (2005) 125119.
  * [25] A. Sekiyama, H. Fujiwara, S. Imada, S. Suga, H. Eisaki, S. I. Uchida, K. Takegahara, H. Harima, Y. Saitoh, I. A. Nekrasov, G. Keller, D. E. Kondakov, A. V. Kozhevnikov, T. Pruschke, K. Held, D. Vollhardt, and V. I. Anisimov: Phys. Rev. Lett. 93 (2004) 156402.
  * [26] S. K. Mo, J. D. Denlinger, H. D. Kim, J. H. Park, J. W. Allen, A. Sekiyama, A. Yamasaki, K. Kadono, S. Suga, Y. Saitoh, T. Muro, P. Metcalf, G. Keller, K. Held, V. Eyert, V. I. Anisimov, and D. Vollhardt: Phys. Rev. Lett. 90 (2003) 186403.
  * [27] L. V. Pourovskii, B. Amadon, S. Biermann, and A. Georges: Phys. Rev. B 76 (2007) 235101.
  * [28] I. L. M. Locht, Y. O. Kvashnin, D. C. M. Rodrigues, M. Pereiro, A. Bergman, L. Bergqvist, A. I. Lichtenstein, M. I. Katsnelson, A. Delin, A. B. Klautau, B. Johansson, I. Di Marco, and O. Eriksson: Phys. Rev. B 94 (2016) 085137.
  * [29] S. Y. Savrasov, G. Kotliar, and E. Abrahams: Nature 410 (2001) 793.
  * [30] L. V. Pourovskii, M. I. Katsnelson, and A. I. Lichtenstein: Phys. Rev. B 72 (2005) 115106.
  * [31] J. Kolorenc, A. B. Shick, and A. I. Lichtenstein: Phys. Rev. B 92 (2015) 085125.
  * [32] L. V. Pourovskii, G. Kotliar, M. I. Katsnelson, and A. I. Lichtenstein: Phys. Rev. B 75 (2007) 235107.
  * [33] J.-Z. Ma, A. van Roekeghem, P. Richard, Z.-H. Liu, H. Miao, L.-K. Zeng, N. Xu, M. Shi, C. Cao, J.-B. He, G.-F. Chen, Y.-L. Sun, G.-H. Cao, S.-C. Wang, S. Biermann, T. Qian, and H. Ding: Phys. Rev. Lett. 113 (2014) 266407\.
  * [34] S. Backes, H. O. Jeschke, and R. Valentí: Phys. Rev. B 92 (2015) 195128.
  * [35] M. Hirayama, T. Miyake, M. Imada, and S. Biermann: arxiv.org arXiv:1511.03757 (2015).
  * [36] A. van Roekeghem, T. Ayral, J. M. Tomczak, M. Casula, N. Xu, H. Ding, M. Ferrero, O. Parcollet, H. Jiang, and S. Biermann: Phys. Rev. Lett. 113 (2014) 266403.
  * [37] A. van Roekeghem, P. Richard, X. Shi, S. Wu, L. Zeng, B. Saparov, Y. Ohtsubo, T. Qian, A. S. Sefat, S. Biermann, and H. Ding: Phys. Rev. B 93 (2016) 245139.
  * [38] A. van Roekeghem and S. Biermann: EPL (Europhysics Letters) 108 (2014) 57003.
  * [39] V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein: J. Phys. Condensed Matter 9 (1997) 767.
  * [40] S. Biermann, F. Aryasetiawan, and A. Georges: Phys. Rev. Lett. 90 (2003) 086402.
  * [41] P. Sun and G. Kotliar: Phys. Rev. Lett. 92 (2004) 196402.
  * [42] J. M. Tomczak, M. Casula, T. Miyake, and S. Biermann: EPL (Europhysics Letters) 100 (2012) 67001.
  * [43] S. Biermann: Journal of Physics: Condensed Matter 26 (2014) 173202.
  * [44] T. Ayral, P. Werner, and S. Biermann: Physical Review Letters 109 (2012) 226401.
  * [45] T. Ayral, S. Biermann, and P. Werner: Physical Review B 87 (2013) 125149\.
  * [46] P. Hansmann, T. Ayral, L. Vaugier, P. Werner, and S. Biermann: Physical Review Letters 110 (2013) 166401.
  * [47] L. Huang, T. Ayral, S. Biermann, and P. Werner: Phys. Rev. B 90 (2014) 195114.
  * [48] L. Boehnke, F. Nilsson, F. Aryasetiawan, and P. Werner: Physical Review B 94 (2016) 201106.
  * [49] T. Ayral, S. Biermann, P. Werner, and L. V. Boehnke: in preparation (2017).
  * [50] J. M. Tomczak, M. Casula, T. Miyake, and S. Biermann: Phys. Rev. B 90 (2014) 165138.
  * [51] M. Casula, P. Werner, L. Vaugier, F. Aryasetiawan, T. Miyake, A. J. Millis, and S. Biermann: Phys. Rev. Lett. 109 (2012) 126408.
  * [52] A. Seidl, A. Görling, P. Vogl, J. Majewski, and M. Levy: Physical Review B 53 (1996) 3764.
  * [53] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz: Phys. Rev. Lett. 49 (1982) 1691.
  * [54] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti: Rev. Mod. Phys. 78 (2006) 865.
  * [55] B. Amadon, S. Biermann, A. Georges, and F. Aryasetiawan: Phys. Rev. Lett. 96 (2006) 066402.
  * [56] N. Xu, P. Richard, A. van Roekeghem, P. Zhang, H. Miao, W.-L. Zhang, T. Qian, M. Ferrero, A. S. Sefat, S. Biermann, and H. Ding: Phys. Rev. X 3 (2013) 011006.
  * [57] P. Werner, M. Casula, T. Miyake, F. Aryasetiawan, A. J. Millis, and S. Biermann: Nature Physics 8 (2012) 331.
  * [58] R. S. Dhaka, Y. Lee, V. K. Anand, D. C. Johnston, B. N. Harmon, and A. Kaminski: Phys. Rev. B 87 (2013) 214516.
  * [59] J. Mansart, P. Le Fèvre, F. m. c. Bertran, A. Forget, D. Colson, and V. Brouet: Phys. Rev. B 94 (2016) 235147.
  * [60] V. Brouet, F. Rullier-Albenque, M. Marsi, B. Mansart, M. Aichhorn, S. Biermann, J. Faure, L. Perfetti, A. Taleb-Ibrahimi, P. Le Fèvre, F. Bertran, A. Forget, and D. Colson: Phys. Rev. Lett. 105 (2010) 087001.
  * [61] R. Poole, R. Leckey, J. Jenkin, and J. Liesegang: Physical Review B 8 (1973) 1401.
  * [62] F. Himpsel, D. Eastman, E. Koch, and A. Williams: Physical Review B 22 (1980) 4604.
  * [63] M. R. Norman: Physical Review B 29 (1984) 2956.
  * [64] F. Aryasetiawan and O. Gunnarsson: Phys. Rev. B 54 (1996) 17564.
  * [65] M. Oshikiri and F. Aryasetiawan: Phys. Rev. B 60 (1999) 10754.
  * [66] J. P. Perdew and A. Zunger: Phys. Rev. B 23 (1981) 5048.
  * [67] H. Jiang, R. I. Gómez-Abal, X.-Z. Li, C. Meisenbichler, C. Ambrosch-Draxl, and M. Scheffler: Computer Physics Communications 184 (2013) 348.
  * [68] B. Holm and U. von Barth: Phys. Rev. B 57 (1998) 2108.
  * [69] M. van Schilfgaarde, T. Kotani, and S. Faleev: Phys. Rev. Lett. 96 (2006) 226402.


[ ]

[] [1]

[]

[cor1]Corresponding author<EMAIL_ADDRESS>
# A Modified Adaptive Genetic Algorithm for Multi-product Multi-period
Inventory Routing Problem

Meysam Mahjoob<EMAIL_ADDRESS>Department of Industrial Engineering,
University of Tehran, Fooman, Rasht  Seyed Sajjad Fazeli
<EMAIL_ADDRESS>Department of Industrial and System Engineering, Wayne
State university, Detroit, MI  Soodabeh Milanlouei
<EMAIL_ADDRESS>Center for Complex Network Research,
Northeastern University, Boston, MA  Leyla Sadat Tavassoli
<EMAIL_ADDRESS>Department of Industrial Manufacturing and
Systems Engineering, University of Texas at Arlington, Arlington, TX

###### Abstract

Recent developments in urbanization and e-commerce have pushed businesses to
deploy efficient systems to decrease their supply chain cost. Vendor Managed
Inventory (VMI) is one of the most widely used strategies to effectively
manage supply chains with multiple parties. VMI implementation asks for
solving the Inventory Routing Problem (IRP). This study considers a multi-
product multi-period inventory routing problem, including a supplier, set of
customers, and a fleet of heterogeneous vehicles. Due to the complex nature of
the IRP, we developed a Modified Adaptive Genetic Algorithm (MAGA) to solve a
variety of instances efficiently. As a benchmark, we considered the results
obtained by Cplex software and an efficient heuristic from the literature.
Through extensive computational experiments on a set of randomly generated
instances, and using different metrics, we show that our approach distinctly
outperforms the other two methods.

###### keywords:

Inventory Routing Problem Genetic Algorithm Vendor Managed Inventory Supply
Chain Management Adaptive Heuristic

##  1 Introduction

Throughout the past few decades, advancements in information and
communications technology have led to novel and innovative businesses in
supply chain management. Vendor Managed Inventory (VMI) is an example of such
business models, representing an important paradigm in which the vendor
(supplier) has full responsibility for controlling the retailer inventory
using the inventory data that the retailer provides. Unlike traditional
inventory management practices where retailers make their own decisions about
the new orders’ size and time, a vendor monitors the retailer’s inventory and
acts as the decision-maker in the VMI model [1]. In this distribution model,
deciding on the inventory level for both supplier and retailer depends on the
timing and the number of deliveries that a retailer needs, which itself is
influenced by the capacity of the vehicles used for delivery. Simultaneous
decision making is a vital component to achieve cost-effectiveness in these
systems. The main advantage of implementing a VMI system is that the vendor
can increase the service level while reducing the distribution costs by using
vehicles effectively. Retailers, on the other hand, can release the resources
they usually use for inventory management. In VMI implementation, a vendor
should solve an Integrated Routing Problem (IRP) to specify a distribution
plan, minimizing the long-run average distribution and inventory costs
throughout the supply chain, and avoiding shortage [2]. The IRP simultaneously
takes into account both vehicle routing and inventory management in a supply
chain system and tries to maximize the overall performance. In the IRP, the
supply chain system consists of a supplier and a number of retailers who are
distributed geographically with certain demand levels. Within a specified time
window, products are shipped from the supplier to retailers using a several
vehicles. The demand of each customer is met by any vehicle per period. In the
VMI system, the supplier determines each retailer’s replenishment policy and
the vehicle routes to deliver the products while ensuring that no shortage
occurs. The IRP attempts to minimize the overall logistics costs to meet
retailers’ demands in a planning horizon. Under this framework, three main
decisions need to be made: 1) when to fulfill the demand of a retailer; 2) the
amount of products to be delivered to a retailer; 3) the optimal routes for
vehicles. The IRP literature can be investigated through the modeling and
methodology perspectives. From the modeling point of view, the literature can
be classified into different groups based on the length of the planning
horizon(single-period, multi-period, and infinite), type of product (single
and multi), and delivery (non-split and split). One of the earliest studies
was done by [3] where authors formulated the IRP as a single-period problem
considering a plant and customers with uncertain demand. The research work by
[4] considered a single product IRP with transshipment where the products
deteriorate during the time of storage or at warehouses. Authors in [5]
considered a multi-period single-product single-vehicle IRPwith a limited
customers to minimize the inventory and transportation cost. An IRP with split
delivery was studied by [6], where the constraint of the VRP that each
customer should be served by exactly one vehicle was relaxed. From the
solution methodology perspective, the researchers have chosen different
approaches toward solving the IRP problem. Although the IRP is NP-Hard, some
researchers attempt to use exact methods to solve the IRP. For example,
authors in [7] developed a Branch-and-Cut (B&C) algorithm to solve a multi-
product multi-vehicle inventory-routing problem. A three-phase exact approach
is developed by [8] where in the first phase, a branch-and-cut algorithm is
applied on extended IRP to obtain an extended lower-bound on cost of IRP. If
the optimal integer solution could not be obtained in the first phase, the
second phase is activated where a heuristic is applied to provide an upper-
bound on the optimal cost. The upper bound will then be used in the third
phase, where a branch-and-cut method is applied on the IRP. While most of the
exact methods for IRP are based on the B&C method, [9] developed a branch-
price algorithm by incorporating a group of valid inequalities, and a labeling
algorithm to solve column generation subproblems, combined with multiple
acceleration techniques.

Due to the complex nature of the IRP problem, many authors developed heuristic
approaches to solve the IRP. A Variable Neighborhood Search (VNS) heuristic is
developed by [10] to solve a multi-period multi-product IRP. They implemented
three different methods to generate the initial solutions. The authors in [11]
proposed a three-phase hybrid heuristic, including inventory allocation,
warehouse clustering, and routing decisions. A Particle Swarm
Optimization–Differential Evolution (PSO-DE) is presented by [12] to solve
green maritime IRP with time window constraints. They implemented a Genetic
Algorithm (GA), Particle Swarm Optimization–Differential Evolution, and basic
PSO to validate the results. A Genetic Algorithm is developed by [13] to solve
a many-to-one distribution network including a depot, an assembly plant and a
number of suppliers. They suggested two types of modifications to their GA
operations to tackle limitations on solving the problem. Authors in [14]
proposed a metaheuristic algorithm to solve the special version of IRP. The
algorithm follows a hierarchy-based consisting of a master problem for route
adjustment and two subproblems for timing and flow decisions.

In this study, we formulate a multi-product multi-period inventory routing
problem (MMIRP) with a supplier, set of retailers (customers) distributed in a
geographical area, and set of heterogeneous vehicles. We develop a Modified
Adaptive Genetic Algorithm (MAGA) where we modify the crossover and mutation
operations based on the chromosome representation matrix. Besides, to
effectively define the GA parameters, we adopt an adaptive approach. The
performance of the MAGA is tested against the Cplex software and a heuristic
approach from the literature. We considered three different metrics to compare
the results. The remainder of this paper is organized as follows: Section 2
present problem definition, assumptions, and mathematical formulation. Section
3 introduces the solution methodology, where we provide details regarding the
proposed heuristic. Section 4 presents instance characteristics, benchmark
methods, and extensive computational experiments. Finally, Section 5 provides
concluding remarks.

##  2 Model Formulation

###  2.1 Problem Definition

In this work, the MMIRP is defined as follows: a given set of vehicles which
are dispatched from a supplier ($w_{0}$), a set of customers which are
distributed in a geographical area; the goal is to determine the efficient
routes for vehicles to satisfy the customers’ demand while minimizing the
fixed fleet cost, inventory, and transportation cost. The detailed assumptions
regarding the model are provided as follows:

  * •

Different types of products are distributed in the network.

  * •

There are no restrictions on the availability of products for the supplier.

  * •

Vehicles are heterogeneous, meaning that they can have different loading
capacities.

  * •

vehicles can visit multiple customers during their trip.

  * •

The demand of customers is deterministic .

  * •

Shortage is not allowed.

###  2.2 Notation

  * •

Sets

    * –

$I$: Set of customers

    * –

$\bar{I}$: Set of customers and supplier where $\bar{I}=I\cup\\{w_{0}\\}$

    * –

$T$: Set of time periods

    * –

$V$: Set of vehicles

    * –

$P$: Set of products types

  * •

Model parameters

    * –

$c_{i,j}$: Travel cost associated with edge $(i,j)\in\bar{I}$.

    * –

$q^{v}$: Capacity of vehicle $v\in V$

    * –

$f_{t}^{v}$: Fixed cost of vehicle $v\in V$ in period $t\in T$

    * –

$s_{i}$: Storage capacity of customer $i\in I$

    * –

$\alpha^{p}$: Weight associated with product $p\in P$

    * –

$d^{p}_{i,t}$: Demand of customer $i\in I$ for product $p\in P$ in time $t\in
T$

    * –

$h^{p}_{i}$: Inventory cost of per unit product $p$ for customer $i$

  * •

Decision variables

    * –

$x_{i,j,t}^{v}$: 1 if the edge $(i,j)\in\bar{I}$ is traversed by vehicle $v\in
V$ in period $t\in T$, and 0 otherwise;

    * –

$y_{i,j,t}^{v,p}$: Amount of product $p\in P$ carried by vehicle $v\in V$ on
edge $(i,j)\in\bar{I}$ .

    * –

$r_{i,t}^{p}$: Amount of on-hand inventory of product $p\in P$ by customer
$i\in I$

###  2.3 Mathematical Formulation

$\displaystyle\text{Min}\sum_{t\in T}\bigg{(}\sum_{i\in I}\sum_{v\in
V}f_{t}^{v}x^{v}_{w_{0},i,t}+\sum_{(i,j)\in\bar{I}}\sum_{v\in
V}c_{i,j}x^{v}_{i,j,t}+\sum_{i\in I}\sum_{p\in P}h^{p}_{i}r_{i,t}^{p}\bigg{)}$
(1) s.t. $\displaystyle\sum_{j\in\bar{I}}x^{v}_{i,j,t}\leq 1\qquad\forall
i\in\bar{I},t\in T,v\in V,$ (2)
$\displaystyle\sum_{j\in\bar{I}}x^{v}_{i,j,t}-\sum_{k\in\bar{I}}x^{v}_{k,i,t}=0\quad\forall
i\in\bar{I},t\in T,v\in V,$ (3) $\displaystyle\sum_{p\in
P}\alpha^{p}y_{i,j,t}^{v,p}\leq q^{v}x_{i,j,t}^{v}\hskip
2.84544pt\forall(i,j)\in\bar{I},t\in T,v\in V,$ (4)
$\displaystyle\sum_{j\in\bar{I}}\alpha^{p}y_{j,i,t}^{v,p}-\sum_{k\in\bar{I}}\alpha^{p}y_{i,k,t}^{v,p}\geq
0\quad\forall i\in\bar{I},t\in T,v\in V,p\in P$ (5) $\displaystyle
r_{i,t-1}^{p}-r_{i,t}^{p}+\sum_{v\in
V}\bigg{(}\sum_{j\in\bar{I}}y_{j,i,t}^{v,p}-\sum_{k\in\bar{I}}y_{i,k,t}^{v,p}\bigg{)}=d^{p}_{i,t}$
$\displaystyle\hskip 99.58464pt\quad\forall i\in\bar{I},t\in T,p\in P$ (6)
$\displaystyle\sum_{p\in P}\alpha^{p}r_{i,t}^{p}\leq s_{i}\qquad\forall i\in
I,t\in T,$ (7) $\displaystyle
x^{v}_{i,j,t}\in\\{0,1\\}\qquad\forall(i,j)\in\bar{I},t\in T,v\in V,$ (8)
$\displaystyle r_{i,t}^{p},y_{i,k,t}^{v,p}\geq 0\qquad\forall i\in\bar{I},p\in
P,t\in T,v\in V$ (9)

The objective function (1) minimizes the fleet fixed, transportation, and
inventory cost. Constraints (2) ensure that each customer is visited at most
once in every period. Connectivity of routes for every vehicle is preserved by
constraints 3. Constraints (4) ensure that the amount of product carried by
each vehicle does not violate the vehicle’s capacity. The sub-tour in any
feasible solution is eliminated by constraint (5). Constraints (6) guarantee
the balance between the demand and on-hand inventory. Also, they ensure that
the demand of customers is satisfied in each period The customers’ storage
capacity violation is prohibited by constraints (7). Finally, constraints (8)
and (9) define the restrictions for the variables.

##  3 Methodology and Algorithm Development

The MMIRP could be solved by commercial solvers. However, to solve the mid and
large-scale instances we need to develop algorithms to efficiently solve the
MMIRP within a reasonable time limit. In the next section we comprehensively
describe the MAGA that we developed to solve the MMIRP.

###  3.1 Modified Genetic Algorithm

The genetic algorithm is a stochastic optimization technique inspired by the
process of natural selection, which is widely applied to solve different
classes of NP-Hard problems [15, 16, 17, 18]. GA maintains a population of
candidate solutions through the selective procedure. GA is initialized by a
set of solutions called population. Each solution in the population is called
a chromosome. Chromosomes progress through successive iterations called
generations. Throughout each iteration, the chromosomes are evaluated by a
fitness function. The fitter chromosomes have higher chances of being selected
for GA operations such as mutation and crossover. The GA operations choose
some parents and produce several offsprings as the new solutions. The new
solutions are then accepted or rejected based on their fitness values as well
as the solutions from the previous iterations to keep the population size
fixed. The GA may converge to the best solution after a certain number of
iterations. To solve the MMIRP using the GA algorithm, we have three major
challenges: 1) The delivery schedule to meet customers’ demand; 2) finding the
best route for each vehicle; and 3) maintaining feasibility such that the
vehicles and customer’s capacity are not violated. Each of these is complex
and difficult to solve, hence, a naive GA may not perform well for this
problem. Therefore, to overcome the complexity, we propose a modified adaptive
GA to efficiently solve the MMIRP.

####  3.1.1 Solution Representation

To represent the solution of the MMIRP, we define a chromosome in the form of
a binary matrix of size ($|I|\times|T|$) where $|I|$ and $|T|$ denote the
number of customers and periods, respectively. The elements $a_{i,t}=1$
indicates that the customer $i$ in period $t$ should be served. The
representation matrix enables us to determine the amount of product for
delivery to each customer in each period. To better demonstrate the advantage
of the representation matrix, we provide an example in the next section.

####  3.1.2 Explanatory Example

We consider a network with a supplier, four customers, two types of product,
four periods, and two vehicles with capacities of 300 and 400 units. Also, the
fixed cost for each of the vehicles is equals to 10 units. The inventory costs
of the two products for each customer in the distribution network are provided
by mini-tables in Figure 1. Also, the transportation cost for each route is
shown on the edges. A feasible solution for this example is illustrated in
Figure 2. Table (a) in this figure is the designed chromosome represented in
the form of a binary matrix. Based on this matrix, customer 1 is served in
periods 1, 3, and 4. Therefore, there is no delivery in period 2, and customer
2’s demand is satisfied by the delivery in period 1. Similarly, the delivery
schedule for the rest of the customers is shown in table (a) in Figure 2.
Tables (b) and (c) in Figure 2 indicate the shipment quantity based on the
weight of products 1 and 2, respectively. For example, the shipment quantity
of products 1 and 2 in period 1 is $1\times(22+3)$ and $2\times(6+26)$
respectively, where 1 and 2 are products’ weight coefficients. Finally, table
(d) 2 shows the total shipment quantity. Once the shipment quantity for
vehicles are calculated, we solve the single-vehicle routing problem by
applying the Lin-Kernighan-Helsgaun heuristic (LKH) [19]. Figure 3 shows the
optimal or near-optimal routes obtained by LKH.

Figure 1: Explanatory example network with a supplier and four customers
Figure 2: A feasible solution for delivery schedules and shipment quantities
represented in a matrix form Figure 3: Optimal or near-optimal routes for
vehicles obtained by LKH method

As it is stated before, there are three types of cost in the distribution
network: 1) fleet’s fixed cost, 2) transportation cost, and 3) inventory cost.
The first two costs are calculated based on the solution obtained in Figure 3.
The inventory cost calculation depends on the period in which the products are
delivered. For example, for the third customer, since the demand for periods 3
and 4 are delivered in period 2, we calculate the inventory cost as
$1\times(27+39)+2\times(12+12)=114$ where the coefficients 1 and 2 indicate
the number of periods that products are stored as inventory.

The feasibility of the chromosome is preserved considering the following
conditions:

  * •

Since the shortage is not allowed, the representation matrix’s elements in the
first column should be 1.

  * •

Split delivery should be avoided.

  * •

In each period, the shipment product should not exceed vehicles’ capacity.

  * •

The amount of delivery to each customer minus the customer’s demand at each
period should not exceed the storage capacity. For instance, suppose that the
storage capacity for customer 4 is violated in period 1, then we schedule a
new delivery one period before the scheduled delivery ($a_{4,3}=1$) to avoid
storage capacity violation.

So, every chromosome is checked based on the above conditions during the GA
process, and the infeasible ones are corrected or dismissed from the solution
pool.

####  3.1.3 Selection

GA’s convergence could be significantly affected by chromosome selection. The
roulette wheel selection was first introduced by [20] to select the
chromosomes for GA operations. Each section of the roulette wheel is assigned
to a chromosome based on the magnitude of its fitness value. The fitness value
of each chromosome is equal to the summation of fleet cost, transportation and
inventory cost. The fitness values of the chromosomes determine their chance
of being selected. We summarize the selection procedure in Algorithm 1. In the
algorithm, the selection and cumulative probabilities of a chromosome $c$ are
referred as $p_{c}$ and $g_{c}$, respectively, and the population size as
$psize$.

Algorithm 1 : Selection

Step 1: Calculate the total fitness F:

$F=\sum_{c=1}^{psize}f_{c}$

Step 2: Calculate the selection probability for each

chromosome:

$p_{c}=\frac{F-f_{c}}{F\times(psize-1)}$

Step 3: Calculate the cumulative probability for each

chromosome:  $g_{c}=\sum_{i=1}^{c}p_{j}\quad c=1,2,...,psize.$

Step 4: Select the chromosome $c$ if :

$g_{c-1}<r\leq g_{c}$

where $r$ is a random number in range (0,1].

####  3.1.4 Crossover

The crossover operator is an essential part of the genetic algorithm, where
two chromosomes (parents) contribute characteristics in creating a new
chromosome (offspring), by randomly exchanging information. Crossover operator
stochastically generates new solutions from an existing population. Crossover
operation increases GA’s ability in searching because it helps inherit and
blend the good properties of the parents in the offspring from the elite
solutions among the population. Since the designed chromosome in this work is
constructed in a matrix form, the crossover operators can be applied on the
rows or columns of the matrix. At each iteration, we randomly select a row or
column from each parent and exchange the selected row or column among the
parents. Figure (4) illustrates the crossover operator on the matrix
chromosome.

Figure 4: The crossover operation on matrix chromosome

####  3.1.5 Mutation

The mutation is another genetic operator that is used to explore new solutions
in the solution space. In some newly formed offspring, some of their genes can
be mutated with a low random probability and they may find characteristics
that do not belong to any of their parents. If the mutation probability is too
high, the searching process will transform into a primitive random search. The
mutation occurs to maintain the diversity within the solution population and
avoid trapping in local optima. The mutation operator designed for this
problem is as follows: Once the roulette wheel algorithm selects a parent, we
randomly select two rows or columns and flip their position. The new
chromosome is then checked by the feasibility conditions and then added to the
pool. Figure 5 shows the mutation procedure.

Figure 5: The mutation operation on matrix chromosome

####  3.1.6 Adaptive GA Operators (AGO)

The performance of GA highly depends on its parameters. In our study, we use
the idea developed by [21]. In their approach, they used the parents’ and
offsprings’ fitness values in each generation to create adaptive crossover and
mutation. The mutation and crossover rates are adjusted based on the returned
fitness value. The better the fitness value returned by any operator, the
higher chance of operator rate. Therefore, the AGO strengthens the well-
performing operators to generate more offsprings. We summarize the AGO steps
in 2. In this algorithm we denote parent size, offspring size, mutation rate,
and crossover rate at iteration $k$ as $pr^{k}_{size}$, $of^{k}_{size}$,
$mr^{k}$, $cr^{k}$. Also, $\bar{F}^{k}_{{pr}_{size}}$ and
$\bar{F}^{k}_{{of}_{size}}$ indicate the average fitness value of parents and
offsprings. We define a stopping criterion as $k_{max}$ which denotes the
number of non-improving steps.

Algorithm 2 : AGO

Initialization:

Set: $k$ $\leftarrow$ $1$,

while $k\leq k_{max}$ :

if $\frac{\bar{F}^{k}_{{pr}_{size}}}{\bar{F}^{k}_{{of}_{size}}}-1\geq 0.1$:

$mr^{k+1}\leftarrow cr^{k}+0.005$, $cr^{k+1}\leftarrow cr^{k}+0.05$

elif $\frac{\bar{F}^{k}_{{pr}_{size}}}{\bar{F}^{k}_{{of}_{size}}}-1\leq 0.1$:

$mr^{k+1}\leftarrow cr^{k}-0.005$, $cr^{k+1}\leftarrow cr^{k}-0.05$

else:

$mr^{k+1}\leftarrow cr^{k}$, $cr^{k+1}\leftarrow cr^{k}$

##  4 Experiments and Results

In this section, we compare the computational performance of the proposed GA
and Constructive Heuristic for Multi-product Inventory Routing Problem
(CHMPIRP) proposed by [22]. The authors employed three performance measures to
evaluate the quality of their heuristic as follows.
Difficulty: is calculated as follows:

$\text{Difficulty (\%)}=\frac{\text{Upper bound}-\text{Lower
bound}}{\text{Upper bound}}$

The lower difficulty percentage shows the tighter bound found by the solver
within the stipulated time-limit. This measure is defined to better represent
the other two measures.
Closeness: The closeness represents the distance between the heuristic’s
objective value (HBV) and the best lower bound obtained by Cplex in
percentage:

$\text{Closeness (\%)}=\frac{\text{HBV}-\text{Lower bound}}{\text{HBV}}$

Large values of Closeness indicate the weak lower bound or poor performance of
heuristic.
Saving: The Saving indicates the distance between the best upper bound
obtained by heuristic’s objective value (HBV) and lower bound of Cplex in
percentage:

$\text{Saving (\%)}=\frac{\text{Upper bound}-\text{HBV}}{\text{Upper bound}}$

Large values of Saving indicates the higher quality of heuristic solution
comparing to the Cplex solver.

They also provided the lower and upper bounds obtained by solving the model
using Cplex 11.0 software package. In our study, all the experiments were
implemented in MATLAB 2019a using a computer with an Intel ®, CPU Core(TM)
i7-9750H, 2.60 GHz, and 8GB RAM. In total 96 instances were generated in a
square grid of size [20, 20], and the supplier at (10,10). The other details
regarding instances are provided in Table 1.

Table 1: Characteristics of the generated instances Parameters | Number of product types
---|---
$|P|=2$ | $|P|=5$
$|I|$ | [5,10,20,30] | [5,10,20,30]
$|T|$ | [5,7] | [5,7]
$|V|$ | [3,5] | [3,5]
$q^{v}$ | $|N|\times 50$ | $|N|\times 100$
$s^{i}$ | 300 | 500
$\alpha^{p}$ | [1,2] | [0.25,0.75,1,1.5,2.5]

###  4.1 Results

Tables (2) and (3) compare the CHMPIRP with the MAGA. The first five columns
indicate the instance number and the size of the instances in terms of the
number of customers, periods, and vehicles. Columns 6 and 7 show the lower and
upper bounds obtained by Cplex. The next two columns represent the best
objective value returned by MAGA and CHMPIRP. Column 10 specifies the
difficulty measure calculated by the lower and upper bounds obtained by Cplex.
The closeness and saving measures calculated for MAGA and CHMPIRP are reported
in the next four columns. The last two columns represent run time in seconds
for MAGA and CHMPIRP, respectively.

Table 2: Computational results for the generated instances Table 3:
Computational results for the generated instances (continue)

In total, 96 instances were solved. In 16 instances the Cplex could not return
any lower or upper bound within the stipulated time limit (2 hours). Also, in
four instances, the Cplex could only return the lower bound. To better compare
the MAGA’s performance with CHMPRIP, we rank the problems based on the
difficulty measure. Figure 6 illustrates the objective values returned by the
two approaches ranked by difficulty percentage. Although the CHMPRIP offers a
better performance in terms of run time, the results show that in almost all
the instances, the MAGA outperformed the CHMPRIP by a significant margin in
terms of the quality of solutions. Figures 7 and 8 represent the performance
of the two algorithms in terms of closeness and saving measures. The MAGA
algorithm was able to provide better closeness in 95 instances by considerable
difference comparing to the CHMPRIP. Also, in 94 instances, the MAGA indicated
superior performance in terms of saving measure compared to the CHMPRIP.

Figure 6: Comparison between the MAGA and CHMPRIP in terms of the objective
functions ranked by difficulty metric  Figure 7: Comparison between the MAGA
and CHMPRIP in terms of closeness metric ranked by difficulty  Figure 8:
Comparison between the MAGA and CHMPRIP in terms of saving metric ranked by
difficulty

###  4.2 Statistical Analysis

To statistically compare the proposed method with the MAGA method, we use the
paired t-test where the null hypothesis is that the mean difference between
two sets of observations is zero. To do so, we solve each of the 96 problems
five times and calculate the average for each problem separately for each
method. We use a paired t-test to compare the means in these two populations,
where observations in one method can be paired with observations in the other
method. The obtained t-statistic is equal to -9.52, resulting in a p-value of
$1.7\text{\times}{10}^{-15}$ which is smaller than the threshold of 0.05,
indicating that there is strong evidence that the mean difference between the
two methods is not equal to zero. With an average difference of -49.47 (95%
CI: -59.78, -39.15), the proposed method generally yields smaller objective
values compared with the MAGA method.

##  5 Conclusion

In this research, we considered a multi-product multi-period inventory routing
problem with a heterogeneous fleet of vehicles. We proposed a modified
adaptive genetic algorithm to efficiently solve a variety of instances. The
approach incorporated different techniques within the genetic algorithm
framework to enhance the performance of the algorithm. Numerical studies were
performed on randomly generated data sets from literature. The efficacy of the
proposed method is benchmarked by comparing it to the commercial solver Cplex
as well as a heuristic from the literature (CHMPRIP). The quality of the
solutions is justified based on three different metrics. Other than a few
instances, the results indicate that MAGA has outperformed the CHMPRIP in all
metrics by a significant margin. Besides, we performed a statistical analysis
to show the significant difference between the two heuristics in terms of
solution quality. Future work could include incorporating green approaches by
either replacing the vehicles with Electric vehicles or incorporating
$\text{CO}_{2}$ emissions into the objective functions. Another extension
could be incorporating uncertainties into the customers’ demand. From an
algorithmic perspective, the use of decomposition algorithms like the column
generation techniques can also be evaluated.

## References

  * [1] Lei Qin, Lixin Miao, Qingfang Ruan, and Ying Zhang.  A local search method for periodic inventory routing problem.  Expert Systems with Applications, 41(2):765–778, 2014.
  * [2] Jean-François Cordeau, Demetrio Laganà, Roberto Musmanno, and Francesca Vocaturo.  A decomposition-based heuristic for the multiple-product inventory-routing problem.  Computers & Operations Research, 55:153–166, 2015.
  * [3] Awi Federgruen and Paul Zipkin.  A combined vehicle routing and inventory allocation problem.  Operations research, 32(5):1019–1037, 1984.
  * [4] A Azadeh, S Elahi, M Hosseinabadi Farahani, and B Nasirian.  A genetic algorithm-taguchi based approach to inventory routing problem of a single perishable product with transshipment.  Computers & Industrial Engineering, 104:124–133, 2017.
  * [5] Claudia Archetti, Luca Bertazzi, Alain Hertz, and M Grazia Speranza.  A hybrid heuristic for an inventory routing problem.  INFORMS Journal on Computing, 24(1):101–116, 2012.
  * [6] Yugang Yu, Haoxun Chen, and Feng Chu.  A new model and hybrid approach for large scale inventory routing problems.  European Journal of Operational Research, 189(3):1022–1040, 2008\.
  * [7] Leandro C Coelho and Gilbert Laporte.  A branch-and-cut algorithm for the multi-product multi-vehicle inventory-routing problem.  International Journal of Production Research, 51(23-24):7156–7169, 2013.
  * [8] Luca Bertazzi, Demetrio Laganà, Jeffrey W Ohlmann, and Rosario Paradiso.  An exact approach for cyclic inbound inventory routing in a level production system.  European Journal of Operational Research, 283(3):915–928, 2020\.
  * [9] Guy Desaulniers, Jørgen G Rakke, and Leandro C Coelho.  A branch-price-and-cut algorithm for the inventory-routing problem.  Transportation Science, 50(3):1060–1076, 2016.
  * [10] Dražen Popović, Milorad Vidović, and Gordana Radivojević.  Variable neighborhood search heuristic for the inventory routing problem in fuel delivery.  Expert Systems with Applications, 39(18):13390–13398, 2012.
  * [11] N Ramkumar, P Subramanian, TT Narendran, and K Ganesh.  A hybrid heuristic for inventory routing problem.  International Journal of Electronic Transport, 1(1):45–63, 2011\.
  * [12] Arijit De, Sri Krishna Kumar, Angappa Gunasekaran, and Manoj Kumar Tiwari.  Sustainable maritime inventory routing problem with time window constraints.  Engineering Applications of Artificial Intelligence, 61:77–95, 2017\.
  * [13] Noor Hasnah Moin, Said Salhi, and NAB Aziz.  An efficient hybrid genetic algorithm for the multi-product multi-period inventory routing problem.  International Journal of Production Economics, 133(1):334–343, 2011\.
  * [14] Zhouxing Su, Zhipeng Lü, Zhuo Wang, Yanmin Qi, and Una Benlic.  A matheuristic algorithm for the inventory routing problem.  Transportation Science, 54(2):330–354, 2020.
  * [15] William Ho, George TS Ho, Ping Ji, and Henry CW Lau.  A hybrid genetic algorithm for the multi-depot vehicle routing problem.  Eng. Appl. Artif. Intell., 21(4):548–557, 2008.
  * [16] Meysam Mahjoob, Seyed Sajjad Fazeli, Soodabeh Milanlouei, Ali Kamali Mohammadzadeh, and Leyla Sadat Tavassoli.  Green supply chain network design with emphasis on inventory decisions.  arXiv preprint arXiv:2104.05924, 2021.
  * [17] Madjid Tavana, Kaveh Khalili-Damghani, Debora Di Caprio, and Zeynab Oveisi.  An evolutionary computation approach to solving repairable multi-state multi-objective redundancy allocation problems.  Neural Computing and Applications, 30(1):127–139, 2018.
  * [18] Leyla Sadat Tavassoli, Nahal Sakhavand, and Seyed Sajjad Fazeli.  Integrated preventive maintenance scheduling model with redundancy for cutting tools on a single machine.  Engineering, Technology & Applied Science Research, 10(6):6542–6548, 2020.
  * [19] Keld Helsgaun.  An effective implementation of the lin–kernighan traveling salesman heuristic.  Eur. J. Oper. Res., 126(1):106–130, 2000.
  * [20] Lawrence Davis.  Applying adaptive algorithms to epistatic domains.  In IJCAI, volume 85, pages 162–164, 1985.
  * [21] KL Mak, YS Wong, and XX Wang.  An adaptive genetic algorithm for manufacturing cell formation.  The International Journal of Advanced Manufacturing Technology, 16(7):491–497, 2000.
  * [22] Nooraddin Dabiri, Mohammad Jafar Tarokh, and Mostafa Setak.  A constructive heuristicfor amulti-productinventoryrouting problem.  Transportation Research, 2(1):11, 2012.