# Non-coherent evolution of closed weakly interacting system leads to
equidistribution of probabilities of microstates

A. P. Meilakhs<EMAIL_ADDRESS>1Centro Atomico Constituyentes, CNEA, Av.
Gral. Paz 1499, San Martin, Buenos Aires, 1650, Argentina

###### Abstract

We introduce a concept of non-coherent evolution of macroscopic quantum
systems. We show that for weakly interacting systems such evolution is
described by a symmetric continuous-time Markov chain. The non-diagonal
elements of its transition rate matrix can be found by Fermi’s golden rule.
Such evolution is time-irreversible and leads to the equidistribution of
probabilities across every state of the system. Further, we search for time
dependence of mean values of occupation number of states and find that under
made assumptions it is governed by Boltzmann collision integral. In this
theory, the non-coherence is the mechanism that transforms time-reversible
unitary evolution into time-irreversible stochastic evolution. Thus we present
the possible solution for the famous time-arrow problem for weakly interacting
systems.

The arrow of time problem, as it is presented in textbooks on statistical
mechanics [1, 2] or specialized books on the nature of time [3, 4], stands as
one of the most fundamental challenges in modern physics. The issue arises
from the time asymmetry in the equations governing the behavior of large
ensembles of particles, in contrast to the symmetric nature of equations
describing individual particles. While discussions regarding the arrow of time
problem typically center around the second law of thermodynamics and entropy
increase, we are interested in exploring other time-asymmetric properties of
large ensembles. Firstly, we want to prove the existence of an equilibrium
state, to which the isolated system strives. Precisely, it is the Gibbs
microcanonical ensemble or the equidistribution of probabilities of
microstates of a macroscopic system. We want to derive the equidistribution of
probabilities in equilibrium from quantum mechanics. Secondly, we seek to
derive the equation that describes the time evolution of gases, namely the
Boltzmann equation. More precisely, we want to derive its irreversible part,
the collision integral.

Probably the most famous approach to the derivation of the Boltzmann equation
from classical mechanics is the BBGKY hierarchy. However, it postulates
molecular chaos as a model for the 2-particle distribution function without
fully elucidating its origins. Another interesting attempt to deduce Boltzmann
equations is made in Refs. [5, 6], where they are treated as mean-field
theory, however only the diffusive, reversible part can be derived this way.
In Ref. [7], authors take the reverse route by not deriving the concept of
probability equidistribution but rather arguing that statistical physics can
be formulated without this assumption.

Other approaches to explain asymmetry postulate low entropy initial states of
the universe [8, 9], which is not completely satisfactory since it is an ad
hoc hypothesis [10]. In contrast, Ref. [11] proposes that physical laws
themselves must be modified to incorporate irreversibility.

Modern approaches propose that irreversibility aemerges from a time-asymmetric
mechanism inherent in quantum mechanics [12]. Typically, decoherence is
suggested as this mechanism [13]. Decoherence theory explains the measurement
in quantum theory as an interaction of a quantum system with media that causes
the loss of coherence in the system [14]. Additionally, non-Markovian quantum
evolution attracted a lot of attention in recent years [15, 16, 17]. It also
focuses on cases where the system dynamics is influenced by the external
world.

We aim to approach the concept of non-coherence from a different perspective.
We observe that it is the coherent states of matter that are exotic, while
non-coherent states occur naturally. Coherent states, such as superconductive
and superfluid states, arise only in specific substances under low
temperatures. Coherent light is generated by lasers. In contrast, natural
light is typically non-coherent, and laser emission gradually loses its
coherence over time without any special conditions. From this, we infer that
non-coherence represents the inherent, natural state of matter.

Rather than focusing on the loss of coherence of a quantum system due to
interaction with its environment, our investigation centers on the environment
itself. We treat it as inherently non-coherent. We aim to comprehend how a
system described by quantum theory but lacking coherence evolves. The
mathematical procedure employed to derive equations for non-coherent light
involves averaging arguments of complex amplitudes of waves. Our objective is
to extend this procedure to encompass general quantum systems and derive the
equation governing non-coherent evolution. We will demonstrate that from the
equation derived in this manner, both the equidistribution of probabilities
and the Boltzmann equation can be readily derived.

The concept of non-coherent transformation was introduced in Ref. [18] to
describe the transition of quantum particles through the interface. It was
introduced to explain the emergence of irreversibility in this case.
Mathematically, a non-coherent transformation is as follows: Let us consider a
set of amplitudes $A_{\alpha}$. The square of the absolute value of each
amplitude represents the probability that the system is in state $\alpha$:
$P_{\alpha}=|A_{\alpha}|^{2}$.At a subsequent moment, the system undergoes
changes, leading to a new set of amplitudes $A^{\prime}_{\alpha}$, which also
represent probability amplitudes. The linear transformation between states
should be unitary to conserve the total probability. We denote matrices and
vectors with specific font: vectors of amplitudes as $\mathcal{A}$,
$\mathcal{A}^{\prime}$, and the unitary transformation matrix as
$\mathcal{U}$. We have

$\mathcal{A}^{\prime}=\mathcal{U}\mathcal{A}.$ (1)

For the amplitude of state $\alpha$ we have

$A^{\prime}_{\alpha}=\sum_{\beta}U_{\alpha\beta}A_{\beta}$ (2)

For probability of the system to be in state $\alpha^{\prime}$, expressed
through amplitudes of the previous state, we have

$P^{\prime}_{\alpha}=|A^{\prime}_{\alpha}|^{2}=\sum_{\beta}|U_{\alpha\beta}|^{2}|A_{\beta}|^{2}+\sum_{\beta\neq\gamma}U^{*}_{\alpha\beta}U_{\alpha\gamma}A^{*}_{\beta}A_{\gamma}$
(3)

We impose the non-coherence condition on the transformation, signifying that
we lack knowledge of the arguments $\varphi$ of the amplitudes in Eq. (1);
only their absolute values are known. To ascertain how the probabilities
transform, we average over all possible values of the arguments:

$P^{\prime}_{\alpha}=\frac{1}{(2\pi)^{N}}\idotsint\limits_{0}^{2\pi}\prod
d\varphi_{\beta}\left[\sum_{\beta}|U_{\alpha\beta}|^{2}|A_{\beta}|^{2}+\sum_{\beta\neq\gamma}U^{*}_{\alpha\beta}U_{\alpha\gamma}A^{*}_{\beta}A_{\gamma}\right]$
(4)

where $N$ is the total number of possible states of the system. We see that
the integral of the sum over $\beta\neq\gamma$ equals zero, and the integral
over the first term equals $(2\pi)^{N}$, which cancels with the denominator.
We introduce $T_{\alpha\beta}=|U_{\alpha\beta}|^{2}$ and write

$P^{\prime}_{\alpha}=\sum_{\beta}T_{\alpha\beta}P_{\beta}.$ (5)

In matrix form it would be

$\mathcal{P}^{\prime}=\mathcal{T}\mathcal{P}.$ (6)

The matrix $\mathcal{T}$ can be expressed as the Hadamard product (elementwise
product) of the matrix $\mathcal{U}$ with its complex conjugate:
$\mathcal{T}=\mathcal{U}^{*}\circ\mathcal{U}$. It follows from the properties
of unitary matrices that the matrix $\mathcal{T}$ is a bistochastic matrix,
meaning its elements are non-negative real numbers, and the sum of the
elements in each row and column equals one. As shown in equation (5), the
values of elements $T_{\alpha\beta}$ can be interpreted as the probability of
transmission from state $\alpha$ to state $\beta$. Bistochastic matrices lack
inverses within the same class of matrices, except for permutation matrices,
which consist solely of ones and zeroes. The physical consequence of this is
that the transformations they describe are irreversible. In Ref. [18], an
experiment is proposed to demonstrate the irreversibility of non-coherent
transmission of light through the interface, thereby validating the presented
theory.

In the discrete case, we observe that the non-coherence condition transforms a
reversible unitary transformation into an irreversible stochastic
transformation. In the continuous case, unitary evolution is described by the
exponential of a skew-Hermitian matrix. In the context of quantum theory, we
typically express this as $\mathcal{U}(t)=\exp(i\mathcal{H}t)$, where
$\mathcal{H}$ is a Hermitian matrix, particularly representing the Hamiltonian
of the system. The continuous analog of a bistochastic transformation is a
symmetric continuous-time Markov chain (CTMC). The time dependence of
transition probabilities is given by $\mathcal{T}(t)=\exp(\mathcal{Q}t)$,
where $\mathcal{Q}$ is the transition rate matrix (Q-matrix) [19]. Its non-
diagonal elements are positive numbers, while the diagonal elements are
negative such that the sums of elements in each column are zero111In
probability texts, the state of the system is often represented by a row, and
the transformation matrix is placed on the right of it. We use the convention
adopted in physical texts, which is to represent the state as a column and
place the transformation matrix to the left of it. This convention results in
a difference in definition: in probability textbooks, the sums of elements in
each row of a Q-matrix are zero.. A general Q-matrix serves as the continuous
analog of a stochastic transformation, while the analog of a bistochastic
transformation is a Q-matrix symmetric over the main diagonal, i.e.,
$Q_{\alpha\beta}=Q_{\beta\alpha}$.

We are going to demonstrate that just as the unitary transformation
universally transitions into a bistochastic transformation under the condition
of non-coherence, unitary evolution transforms into symmetric stochastic
evolution, albeit exclusively in the case of a closed macroscopic system. Both
conditions are crucial: the closed nature of the system ensures a discrete
spectrum, and we focus solely on interactions within the system. We consider a
macroscopic system where a continuous spectrum with a well-defined density of
states serves as a suitable approximation. The CTMC turns out to be as good an
approximation of the evolution of the system as the density of states is a
good approximation of its spectrum. We also imply the condition of weak
interaction since strongly interacting systems are beyond the applicability of
the methods that we use.

We start similarly to classic time-dependent perturbation theory [20]. We seek
for a wave function of a many-body system in the form

$\Psi=\sum_{\alpha}A_{\alpha}|\alpha\rangle,$ (7)

where $|\alpha\rangle$ are states of non-interacting system, and $A_{\alpha}$
are probability amplitudes of system to be in state $|\alpha\rangle$.
Standartly we get the following equation for $A_{\alpha}$:

$i\hbar\sum_{\alpha}\dot{A}_{\alpha}|\alpha\rangle=\sum_{\beta}\exp\left(\frac{i}{\hbar}(E_{\alpha}-E_{\beta})t\right)\hat{V}A_{\beta}|\beta\rangle.$
(8)

We introduce $V_{\alpha\beta}=\langle\alpha|\hat{V}|\beta\rangle$ and get the
equation for amplitudes

$\dot{A}_{\alpha}=-\frac{i}{\hbar}\sum_{\beta}V_{\alpha\beta}\exp\left(\frac{i}{\hbar}(E_{\alpha}-E_{\beta})t\right)A_{\beta},$
(9)

We note that this equation describes unitary evolution of amplitudes since
$V_{\alpha\beta}e^{i(E_{\alpha}-E_{\beta})t/\hbar}=\left(V_{\beta\alpha}e^{i(E_{\beta}-E_{\alpha})t/\hbar}\right)^{*}$.
Also, interaction occurs between different states so diagonal elements equal
zero $V_{\alpha\alpha}=0$.

Now we are departing from the standard approach. Let us assume amplitudes are
slowly changing function over time compared to fast oscillating exponent
$\exp(i(E_{\alpha}-E_{\beta})t/\hbar)$ for most values of $\alpha$ and
$\beta$. We want to distinguish slow oscillating terms with $E_{\alpha}\approx
E_{\beta}$ that contribute most to the change of amplitudes. To do that we
integrate over some time range $\Delta t$ that is small compared to the
typical timescale of change of $A_{\alpha}$, however
$(E_{\alpha}-E_{\beta})\Delta t/\hbar\gg 1$ for most values of $\alpha$ and
$\beta$. We will show that if such a time scale exists the results we obtain
will not depend on the value of $\Delta t$. On the left side of the equation,
we obtain $A(\Delta t)-A(0)$. On the right side

$\displaystyle\int\limits_{0}^{\Delta
t}\exp\left(\frac{i}{\hbar}(E_{\alpha}-E_{\beta})t\right)A\approx
A\int\limits_{0}^{\Delta
t}\exp\left(\frac{i}{\hbar}(E_{\alpha}-E_{\beta})t\right)=-\frac{i\hbar}{E_{\alpha}-E_{\beta}}\left(e^{\frac{i}{\hbar}(E_{\alpha}-E_{\beta})\Delta
t}-1\right)A$ (10)

We will define the expression before $A$ in the right side of the equation as
$\chi(\Delta E)$, where $\Delta E=|E_{\alpha}-E_{\beta}|$. We divide both
sides by $\Delta t$. On the left we get

$\displaystyle\frac{A(\Delta t)-A(0)}{\Delta t}\approx\dot{A},$ (11)

but here we have a time derivative of smoothed $A$, which we have truncated
the fast oscillating terms from. For the right side, we will consider two
cases:

$\displaystyle\Delta E\Delta t\gg\hbar=>\chi(\Delta E)\approx 0,$
$\displaystyle\Delta E\Delta t\ll\hbar=>\chi(\Delta E)\approx 1.$ (12)

We proceed by exchanging the $\chi(\Delta E)$ to a new function
$\overline{\chi}(\Delta E)$ which equals one in a region of energies near
$\Delta E=0$ and zero outside of this region. Since we use $\chi$ to define
integral transformation of square integrable function, particulary in formula
$\sum_{\beta}\chi(E_{\alpha}-E_{\beta})$, we impose the condition $\int
dE|\chi|^{2}=\int dE|\overline{\chi}|^{2}$ which yields

$\displaystyle\Delta E\Delta t>$
$\displaystyle\pi\hbar=>\overline{\chi}(\Delta E)=0$ $\displaystyle\Delta
E\Delta t<$ $\displaystyle\pi\hbar=>\overline{\chi}(\Delta E)=1.$ (13)

We will discuss later what we have dismissed with this simplification.

We come to the equation for smoothed amplitudes:

$\dot{A}_{\alpha}=-\frac{i}{\hbar}\sum_{\beta}V_{\alpha\beta}A_{\beta},$ (14)

where $E_{\alpha}-E_{\beta}<\Delta E$.

We want to write the equation in the matrix form

$\frac{d}{d\,t}\begin{pmatrix}A_{1}\\\ A_{2}\\\
...\end{pmatrix}=-\frac{i}{\hbar}\begin{pmatrix}V_{11}&V_{12}&...\\\
V_{21}&V_{22}&...\\\ ...&...&\end{pmatrix}\begin{pmatrix}A_{1}\\\ A_{2}\\\
...\end{pmatrix}$ (15)

We will denote vectors rows and matrices by special font, and rewrite the
previous equation as

$\dot{\mathcal{A}}=-\frac{i}{\hbar}\mathcal{V}\mathcal{A}$ (16)

We note that matrix $\mathcal{V}$ is hermitian so equation (16) describes
unitary evolution of the system.

The solution is

$\mathcal{A}(t)=\exp\left(-\frac{i}{\hbar}\mathcal{V}t\right)\mathcal{A}(0)$
(17)

The time-derivative of propability of system being in state $\alpha$ is
$\dot{P}_{\alpha}=\dot{A}^{*}_{\alpha}A_{\alpha}+A^{*}_{\alpha}\dot{A}_{\alpha}$.
We use approximate derivative defined by equation (11). We need to find
$\mathcal{A}(\Delta t)$. Since $\Delta t$ is small compared to timescale of
change of $\mathcal{A}$, so we can expand exponent to the Taylor series with
first-order precision in $\Delta t$:

$\mathcal{A}(\Delta t)=\left(\mathcal{E}-\frac{i}{\hbar}\mathcal{V}\Delta
t\right)\mathcal{A}(0)$ (18)

where $\mathcal{E}$ is an unit matrix. We write it element wise so we can see
first non-vanishing term for each element:

$A_{\alpha}(\Delta
t)=A_{\alpha}(0)-\frac{i}{\hbar}\sum_{\beta}V_{\alpha\beta}\Delta
tA_{\beta}(0)$ (19)

And equivalent expressin we can write for $A^{*}_{\alpha}$. Then for
derivative, to the first non-vanishing term in $\Delta t$

$\dot{A}_{\alpha}(\Delta
t)=-\frac{i}{\hbar}\sum_{\beta}V_{\alpha\beta}A_{\beta}(0)+\frac{1}{\hbar^{2}}\left(\sum_{\beta}|V_{\alpha\beta}|^{2}\right)\Delta
tA_{\alpha}(0),$ (20)

where we have used that $\mathcal{V}$ is Hermitian.

Multiply expressions (19, 20) and imply non-coherence condition (4) that
guarantees that $A^{*}_{\alpha}A_{\beta}=\delta_{\alpha\beta}|A_{\alpha}|^{2}$
to obtain

$\dot{P}_{\alpha}=-\frac{2}{\hbar^{2}}\left(\sum_{\beta}|V_{\alpha\beta}|^{2}\right)\Delta
tP_{\alpha}+\frac{2}{\hbar^{2}}\sum_{\beta}|V_{\alpha\beta}|^{2}\Delta
tP_{\beta}$ (21)

We introduce Q-matrix defining

$\displaystyle\alpha\neq\beta,\ \ \
Q_{\alpha\beta}=\frac{2}{\hbar^{2}}|V_{\alpha\beta}|^{2}\Delta t$
$\displaystyle Q_{\alpha\alpha}=-\frac{2}{\hbar^{2}}\Delta
t\left(\sum_{\beta}|V_{\alpha\beta}|^{2}\right).$ (22)

We can see that it satisfies the conditions for the transition rate matrix.
Non-diagonal elements are positive numbers and diagonal elements are negatives
such that the sums of elements in each column are zero. The evolution of
probabilities is governed by the equation

$\dot{\mathcal{P}}=\mathcal{Q}\mathcal{P}.$ (23)

It is a stochastic evolution that converges to the stationary distribution for
large values of t. The Q-matrix is derived from the Hermitian matrix, ensuring
its symmetry, which aligns with the principle of detailed balance. Symmetry
guarantees that the stationary distribution is an equiprobability of states.
At the limit $t\to\infty$ probability of every state becomes equal [19]. That
means the distribution to which the non-coherent evolution of a closed
macroscopic system converges is the Gibbs microcanonical ensemble.

Let’s address the issue with formula (22): its dependence on the arbitrarily
chosen parameter $\Delta t$. We will now demonstrate that for macroscopic
systems, or as we define them, systems for which the density of states is a
good approximation, the system’s dynamics do not rely on the choice of $\Delta
t$.

We begin by utilizing an approximation of a continuous spectrum and, near a
given state $\beta$, we express quantum numbers as two components. The first
is $\lambda_{T}$ such that little change in its value leaves energy unchanged,
so it is tangential to the energy surface. The second is $\lambda_{N}$ which
change results in change in energy, so we write energy as a function of quntum
numbers in form
$E(\lambda)=E(\lambda_{0})+\lambda_{N}\mathrm{grad}_{\lambda}E$. Also we
recall expression (13) and substitute $\Delta t=\pi\hbar/\Delta E$. We obtain:

$\frac{2}{\hbar^{2}}\sum_{\beta}|V_{\alpha\beta}|^{2}\Delta
tP_{\beta}=\frac{2\pi}{\hbar\Delta
E}\sum_{\lambda_{T},\lambda_{N}}|V_{\alpha\beta}|^{2}P_{\lambda_{T},\lambda_{N}}$
(24)

If $P_{\lambda_{N}}$ is a smooth function of energy at the point of interest,
in a short interval we can approximate it as a constant, so summation over
$\lambda_{N}$ gives the number of states inside the energy region – $N$:

$\frac{2\pi}{\hbar\Delta
E}\sum_{\lambda_{T},\lambda_{N}}|V_{\alpha\beta}|^{2}P_{\lambda_{T}}=\frac{2\pi
N}{\hbar\Delta
E}\sum_{\lambda_{T}}|V_{\alpha\beta}|^{2}P_{\lambda_{T}}=\frac{2\pi}{\hbar}\sum_{\lambda_{T}}|V_{\alpha\beta}|^{2}P_{\lambda_{T}}\nu(E)$
(25)

where $\nu(E)$ is the dencity of states $N/\Delta E$.

We can reason differently. We look at the derived formulae for the transition
rates (21,22). When $\Delta t\to\infty$, the region of energies shrinks to
infinitesimal ($dE$), but the value of the transmission rate grows to infinity
(delta function). For transmission between states $\alpha$ and $\beta$ we get

$Q_{\alpha\beta}=\frac{2\pi}{\hbar}|V_{\alpha\beta}|^{2}\nu(E_{\beta})\delta(E_{\alpha}-E_{\beta}).$
(26)

So we have arrived at the well-known Fermi’s golden rule.

It is a good moment to remember what we have dropped by using approximation
(13). At large $\Delta t$ each term in summation (14) outside of a small
region of energies $\Delta E$ gives a small contribution, however, there are a
lot of such terms so just discarding them all is not accurate. We can write:

$\frac{1}{\Delta
t}\frac{i\hbar}{(E_{\alpha}-E_{\beta})}=\frac{i}{\pi}\frac{\Delta
E}{(E_{\alpha}-E_{\beta})}\underset{\Delta
t\to\infty}{=}\frac{i}{\pi}\frac{dE}{(E_{\alpha}-E_{\beta})}.$ (27)

And we can exchange summation to integration over energy. This leads us to the
expression for the perturbation of energies and states by interaction.
Notably, as we approach the limit $\Delta t\to\infty$, our formulas align with
those of well-established perturbation theory. Therefore, we do not contradict
existing theories but aim to elucidate the physical interpretation of their
mathematical formalism.

We see that with approximation (13) we disregard perturbation of states by
interaction. We can speculate that considering interaction, we will still have
stochastic evolution described by equation (23) but with renormalized energies
of states and matrix elements of interaction. However, we can not prove it and
will proceed to assume interaction is weak so we can disregard perturbation of
states by interaction.

In weak interaction limit states of the system may be expressed by occupation
numbers of states of one-particle Hamiltonian. Let us consider a typical case:
transmissions between levels by fermions with emission or absorption of
bosons. For example, it could be electrons and phonons. The matrix element has
the form $V_{\alpha\beta}\hat{c}^{\dagger}_{\alpha}\hat{c}_{\beta}\hat{a}$,
where $\hat{c},\hat{a}$ are annihilation operators for fermions and bosons
respectively.

Let us first consider just two fermionic states and the transitions between
them. We denote states of a noncoherent system by angular brackets on both
sides, in this notation we consider transitions $\langle
0,1,n-1\rangle\leftrightarrow\langle 1,0,n\rangle$. We use formula (21) with
$\alpha=\langle 0,1,n-1\rangle,\beta=\langle 1,0,n\rangle$ and
$q=2\pi|V|^{2}_{\alpha\beta}/\hbar$:

$\displaystyle\dot{P}_{\langle 0,1,n-1\rangle}=-nqP_{\langle
0,1,n-1\rangle}+nqP_{\langle 1,0,n\rangle}$ $\displaystyle\dot{P}_{\langle
1,0,n\rangle}=-nqP_{\langle 1,0,n\rangle}+nqP_{\langle 0,1,n-1\rangle}.$ (28)

We want to find how the mean occupation number (which is a one-particle
distribution function) of the first state $f_{1}$ varies in time. We write
down the expression $f_{1}$ through probabilities

$\displaystyle f_{1}=$ $\displaystyle P_{\langle 0,1,0\rangle}\times
0+P_{\langle 0,1,1\rangle}\times 0+P_{\langle 0,1,2\rangle}\times 0+...$
$\displaystyle+P_{\langle 1,0,0\rangle}\times 1+$ $\displaystyle P_{\langle
1,0,1\rangle}\times 1+P_{\langle 1,0,2\rangle}\times 1+...$ (29)

Vertically aligned are states between which the transitions occur.

We take time derivative from both sides and using (28) we get

$\displaystyle\dot{f}_{1}=(-qP_{\langle 1,0,1\rangle}+qP_{\langle
0,1,0\rangle})+(-2qP_{\langle 1,0,2\rangle}+2qP_{\langle
0,1,1\rangle})+...=\sum_{n}nq(P_{\langle 0,1,n-1\rangle}-P_{\langle
1,0,n\rangle}).$ (30)

To deduce one-particle distribution function from many-particle probabilities
we need to make a simplifying assumption. We assume that the probabilities of
different occupation numbers of different states are independent of each
other. That means $P_{\langle i,j,n\rangle}=P_{\langle
i,...\rangle}P_{\langle...,j\rangle}P_{\langle n\rangle}$ . We substitute this
expression into (31) and obtain

$\displaystyle\dot{f}_{1}=\sum_{n}nq(P_{\langle
0,...\rangle}P_{\langle...,1\rangle}P_{\langle n-1\rangle}-P_{\langle
1,...\rangle}P_{\langle...,0\rangle}P_{\langle
n\rangle})=q\sum_{n}n\left[(1-f_{1})f_{2}P_{\langle
n-1\rangle}-f_{1}(1-f_{2})P_{\langle n\rangle})\right]=$
$\displaystyle=q\left[(1-f_{1})f_{2}\sum_{k=0}(k+1)P_{\langle
k\rangle}-f_{1}(1-f_{2})\sum_{k=0}kP_{\langle
k\rangle}\right]=q\left[(1-f_{1})f_{2}(N+1)-f_{1}(1-f_{2})N\right].$ (31)

Here $N$ denotes the mean occupation number of phonons. Summation over all
possible transitions yeilds the Boltzmann collision integral. With analogous
arguments we can derive collision integral for other cases, such as three-
phonon processes.

Thus we have deduced that the limiting point of evolution of closed
macroscopic system is the Gibbs microcanonical ensemble which is the basis for
all the equilibrium statistical mechanics. We also derived the collizion part
of the Boltzmann equation which is the fundamental equation of physical
kinetics. This was achieved by postulating the incoherence of the system in
consideration, which practically means averaging arguments of complex
amplitudes for different states of the system. We drew an analogy with non-
coherent light, where averaging over the arguments of waves provides correct
and testable predictions.

In addition, we want to present two mathematical arguments, which do not
rigorously prove the correctness of such averaging, but at least make it
believable. First, the presented theory describes transitions to some energy
regions, which contain a lot of very close states. There is no a-priory reason
why phases of these states should be correlated, so averaging over arguments
is natural. Second, we describe the transition over time $\Delta t$, to a
region of frequencies $\Delta E/\hbar$. We have a condition $\Delta t\Delta
E/\hbar>1$, which precisely aligns with the condition for the loss of
coherence in wave optics. Consequently, our theory applies to time intervals
larger than the coherence time.

We notice that such arguments apply not only to the macroscopic closed system
but also, for example, to a single atom emitting photon into space. The Fermi
rule was originally derived for such a situation. While an atom has a discrete
spectrum the atom together with an electromagnetic field has a continuous
spectrum so arguments hold. If we think the quantum system loses coherence in
interaction with the macroscopic body, in the context of the presented theory
vacuum is a macroscopic body. The difference between an isolated system and an
atom emitting into space is photon departs from the atom and will not be
reabsorbed, so the symmetry of forward and backward transitions is disrupted.

We have basiclly re-derived the Fermi rule with an application to statistical
mechanics. Fermi rule calculates probabilities of transitions and not
amplitudes of transitions so it always was about stochastic evolution rather
than unitary evolution. Emphasizing this aspect of the well-known formula and
associating it with non-coherence, we provide a novel perspective on its
interpretation. Additionally, we have calculated the diagonal elements for the
Q-matrix, which further confirms our understanding of transmissions between
states as the manifestation of stochastic evolution.

Other than the non-coherence condition, an important approach used in the
presented paper is in the derivation of transmission rates. Traditional
quantum perturbation theory methods involve integration over an infinite
period of time. However, in our work, we sought to localize transitions
somewhere on a timeline. To achieve this, we integrated over a finite region
of time $\Delta t$ and later show that in the limit $\Delta t\to\infty$ we
recover the well-known formulae of perturbation theory (26, 27). Having finite
value throughout the derivation is crucial. First, it allows us to use
approximation (11). Then we use a discrete structure of levels of the system
to apply the CTMC theory to our situation.

We can not also consider $\Delta t$ small after equation (14) is derived.
Finiteness of $\Delta t$ allows us to keep the second summand in Eq. (20)
which in turn allows to find the diagonal term in Q-matrix (22).

We have derived the expression (23) where transitions between states are
assumed to be instantaneous. However, our derivation allows us to see that it
is only negligibly small compared to the characteristic timescale of the
evolution of amplitudes. When we take a limit $\Delta t\to\infty$ we also
would like to know, compared to which quantity it is big. After all $\Delta t$
is a dimensional parameter and should be compared to something with a
dimension of time. Hence we want to estimate the transition time $\tau$.

For simplicity let us think about the emission of atom rather than transitions
in macroscopic system. First,we observe that photons with energy differences
smaller than $\Delta E$ are effectively ”glued” together or indistinguishable
in our formalism. Hence, photons with energies less than $\Delta E$ are
practically indistinguishable from an absence of any photons. On a timescale
$\Delta t$ we only can see a photon with energy $\hbar\omega>\hbar/\Delta t$
hence $T<\Delta t$, where $T$ is a wave period. On the other hand, we can not
separate in time events that are separated by less than $\Delta t$. The
transition time is the minimal time interval that still allows the detection
of the emission $\tau=inf\\{\Delta t:\,T<\Delta t\\}$. This means $\tau=T$,
the transition time is equal to the period of the wave emitted, which makes
perfect sense.

Thus we have found the limit of applicability of the presented theory: such a
timescale $\Delta t$ should exist that $\tau\ll\Delta t\ll T_{A}$, where
$T_{A}$ is a characteristic timescale of change of amplitudes.

It looks like we have ”lost” waves of small frequencies $\hbar\omega<\Delta E$
in our consideration. We suggest they contribute to the diffusion term in the
Boltzmann equation, and should be treated separately (see Ref. [5]).

## References

  * [1] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics, Vol. 5: Statistical Physics Part 1 (Pergamon, New York, 1980) [Chapter 1, §8].
  * [2] D. Kondepudi, I. Prigogine Modern Thermodynamics: From Heat Engines to Dissipative Structures (Wiley, New York, 2014) [Chapter 3, §4]
  * [3] H. D. Zeh, The Physical Basis of the Direction of Time (Springer-Verlag, Berlin, 1999).
  * [4] A. L. Kuzemsky, The Mystery Of Time: Asymmetry Of Time And Irreversibility In The Natural Processes (World Scientific, New Jersey, 2022)
  * [5] F. Golse, Commun. Math. Phys. 343, 165–205 (2016).
  * [6] F. Golse, arXiv:2201.02005[math-ph]
  * [7] S. Popescu, A. J. Short, A. Winter Nature Phys. 2, 754–758 (2006).
  * [8] R. M. Wald, Stud. Hist. and Philosophy of Sci. 37, 394-398 (2006).
  * [9] S. M. Carroll and J. Chen, Gen. Relativ. Gravit. 37, 1671 (2005).
  * [10] E. G. D. Cohen and T. H. Berlin, Physica (Amsterdam) 26, 717 (1960).
  * [11] G. C. Ghirardi, A. Rimini, and T. Weber, Phys. Rev. D 34, 470 (1986).
  * [12] L. Maccone, Phys. Rev. Lett. 103, 080401 (2009).
  * [13] L. S. Schulman, Time’s Arrow and Quantum Measurement (Cambridge University Press, Cambridge, England, 1997).
  * [14] E. Joos, H. D. Zeh, C. Kiefer, D. Giulini, J. Kupsch, I.-O. Stamatescu Decoherence and the Appearance of a Classical World in Quantum Theory (Springer, Berlin, 2003).
  * [15] H. Spohn, Rev. Mod. Phys. 52, 569 (1980).
  * [16] M. Merkli, Ann. Phys. 412, 167996 (2020).
  * [17] D. Chruscinski, Phys. Reports 992, 1-85 (2022).
  * [18] A. P. Meilakhs, arXiv:2401.00059 [cond-mat.other]
  * [19] J. R. Norris, Markov Chains (Cambridge University Press, Cambridge, 1997)
  * [20] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics, Vol. 3: Quantum mechanics (Butterworth-Heinemann, Oxford , 1981) [Chapter VI, §40].


AAffil[arabic] ANote[fnsymbol]

# An Algorithm for Adversary Aware Decentralized Networked MARL

Soumajyoti Sarkar This work has been done while the author was at Arizona
State University and
prior to joining Amazon  CaliforniaUnited States<EMAIL_ADDRESS>
###### Abstract.

Decentralized multi-agent reinforcement learning (MARL) algorithms have become
popular in the literature since it allows heterogeneous agents to have their
own reward functions as opposed to canonical multi-agent Markov Decision
Process (MDP) settings which assume common reward functions over all agents.
In this work, we follow the existing work on collaborative MARL where agents
in a connected time varying network can exchange information among each other
in order to reach a consensus. We introduce vulnerabilities in the consensus
updates of existing decentralized MARL algorithms where some agents can
deviate from their usual consensus update, who we term as adversarial agents.
We then proceed to provide an algorithm that allows non-adversarial agents to
reach a consensus in the presence of adversaries under a constrained setting.

##  1\. Introduction

Multi-agent reinforcement learning constitutes a reinforcement learning
scenario where multiple agents participate to jointly learn an optimal policy
and where optimality is defined in terms of some objective (Busoniu et al.,
2008). In a collaborative setting, the agents have to work together with a
goal to optimize a shared reward metric (Kok and Vlassis, 2006). Multi-agent
environments are inherently non-stationary since the other agents are free to
change their behavior as they also learn and adapt and this makes designing
algorithms for MARL more complex that single agent systems. Our model differs
from the traditional collaborative and fully co-operative multi-agent systems
in that the agents in our system are heterogeneous and can have different
local reward functions. We work on sequential decision-making problems in
which the agents repeatedly interact with their environment and try to jointly
optimize the long-term reward they receive from the system, which depends on a
sequence of decisions made by them and the signals shared by other agents in
the network.

Existing value-factorized based MARL approaches perform well in various multi-
agent cooperative environments under the centralized training and
decentralized execution (CTDE) scheme, where all agents in a team are trained
together by the centralized value network and each agent executes its policy
independently (Lyu et al., 2021; Lowe et al., 2017). One of the existing
issues in such systems is that in the centralized training process, the
environment for the team is partially observable and non-stationary. The
observation and action information of all the agents cannot represent the
global states. What follows is that the existing methods perform poorly and
are highly sample inefficient. In these regimes, regret minimization is a
promising approach as it performs well in partially observable and fully
competitive settings. However, it tends to model others as opponents and thus
cannot work well under the CTDE scheme. In our work, we set out to respect two
constraints: first, we want the agents to be fully decentralized even in their
training and second, we follow the model free actor critic type algorithms for
our collaborative MARL setup as opposed to Q-learning type algorithms or
centralized training with model based approaches (Willemsen et al., 2021).

The major difference in the kind of algorithms we follow in this paper with
that of CTDE techniques is that our work follows the studies on consensus
based distributed optimization (Nedic and Ozdaglar, 2009) where the agents
exchange parameters instead of their actions and states with other agents
(Lowe et al., 2017). There is no centralized training and each agent performs
the critic step independently. These consensus based systems have been popular
in many fields like sensor networks, social learning, co-ordination of
vehicles and co-ordinated synchronous distributed optimization to name a few.
However, most of these studies evade the risk of adversaries or agents that
deviate from consensus updates. For example, in studies where humans
participate in decision making and social learning over time (Celis et al.,
2017; Sarkar et al., 2020), presence of adversarial nodes could disrupt the
learning dynamics. The study in (Figura et al., 2021) showed that a single
adversarial agent can persuade all the other agents in the network to
implement policies that optimize an objective that it desires. It then becomes
important to answer whether we can modify the existing consensus based
decentralized MARL algorithms to respond to the presence of malicious agents
who deviate from the consensus rule updates or whether the non-malicious
agents can still converge to the optimal solution in the presence of
adversaries.

We study the setup where the rewards or the data for the agents are not
corrupted but instead where the adversarial agents do not follow the consensus
updates to obtain a joint optimal policy. We consider the scenario where there
could exist more than a single adversarial agent in our setup. Section 2
provides the technical preliminaries of our networked MARL setup and Section 3
discusses the adversarial settings and the adversary aware consensus MARL
algorithm.

##  2\. Decentralized Networked MARL

Our work is heavily influenced by the research on decentralized multi-agent
reinforcement learning done in (Zhang et al., 2018). The authors develop a
consensus based reinforcement learning algorithm in a fully decentralized and
networked setting. We primarily extend their work in a networked setting where
malicious agents are present under a constrained setup. Our networked
multiagent MDP is constructed as follows: we have a state space $\mathcal{S}$
shared by all agents $\mathcal{N}$ in a network, such that $\mathcal{A}^{i}$,
$i\in|\mathcal{N}|$ is the action of agent $i$. Let
$\mathcal{A}=\prod_{i=1}^{N}\mathcal{A}^{i}$ be the joint action space of all
agents. Following this,
$\mathcal{R}^{i}:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is the
local reward function of agent $i$ and
$P:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow[0,1]$ is the state
transition probability of the MDP. We assume that the states are observable
globally in order to ensure co-operation among agents, but the rewards are
observed locally as a function of the neighbors of the agent. At time $t$, an
agent $i$ chooses an action $a^{i}_{t}$ given state $s_{t}$, following its own
local policy $\pi^{i}:\mathcal{S}\times\mathcal{A}^{i}\rightarrow[0,1]$, where
$\pi^{i}(s,a^{i})$ represents the probability of choosing action $a^{i}$ at
state $s$. In this setup, the joint policy of all agents is given by
$\pi(s,a)=\prod_{i\in\mathcal{N}}\pi^{i}(s,a^{i})$.

For an agent $i$, the local policy is parameterized by $\pi^{i}_{\theta}$
where $\theta^{i}\in\Theta^{i}$ are the parameters and the joint policy is
given by $\Theta=\prod_{i=1}^{N}\Theta^{i}$. As is standard in the assumptions
of actor-critic algorithms with function approximation, the policy function
$\pi_{\theta}^{i}(s,a^{i})>0$ for any $\theta^{i}\in\Theta^{i}$. We assume
that $\pi_{\theta^{i}}^{i}(s,a^{i})$ is continuously differentiable with
respect to the parameter $\theta^{i}$ over $\Theta^{i}$. For all
$\theta\in\Theta$, the transition matrix of the Markov chain
$\\{s_{t}\\}_{t\geq 0}$ induced by policy $\pi_{\theta}$, for any
$s,s^{\prime}$ $\in\mathcal{S}$ is given by

(1) $P^{\theta}(s^{\prime}|s)=\sum_{a\in\mathcal{A}}\pi_{\theta}(s,a)\cdot
P(s^{\prime}|s,a)$

It comes along with the assumption that $\\{s_{t}\\}_{t\geq 0}$ is irreducible
and aperiodic under any $\pi_{\theta}$, and the chain of the state-action pair
$\\{(s_{t},a_{t})\\}$ has a stationary distribution
$d_{\theta}(s)\cdot\pi_{\theta}(s,a)$ for all pairs.

Objective of the MARL: The collective objective of the agents is to
collaborate and find a policy $\pi_{\theta}$ such that it maximizes the
average global long term rewards while only utilizing information that is
local to the agents in the network. In that context, let $r_{t+1}^{i}$ denote
the reward received by agent $i$ at time $t$. Then the goal of all agents
collectively is to optimize the following objectives:

(2)
$\begin{split}max_{\theta}J(\theta)&=lim_{T}\frac{1}{T}\mathbb{E}\Big{(}\sum_{t=0}^{T-1}\frac{1}{N}\sum_{i\in\mathcal{N}}r_{t+1}^{i}\Big{)}\\\
&=\sum_{s\in\mathcal{S},a\in\mathcal{A}}d_{\theta}(s)\pi_{\theta}(s,a)\cdot\overline{R}(s,a).\end{split}$

where $\overline{R}(s,a)$ = $\frac{1}{N}\cdot\sum_{i\in\mathcal{N}}R^{i}(s,a)$
is the globally averaged reward function. Following this, we have
$\overline{R}(s,a)=\mathbb{E}[\overline{r}_{t+1}|s_{t}=s,a_{t}=a]$ where
$\overline{r}_{t}=\frac{1}{N}\sum_{i\in\mathcal{N}}t_{t}^{i}$. In that context
of the symbols, the global state value and action value functions can be
denoted by
$\mathcal{Q}_{\theta}(s,a)=\sum_{t}\mathbb{E}[\overline{r}_{t+1}-J(\theta)\ |\
s_{0}=s,a_{0}=a,\pi_{\theta}]$ and
$\mathcal{V}_{\theta}(s)=\sum_{a\in\mathcal{A}}\pi_{\theta}(s,a)\mathcal{Q}_{\theta}(s,a)$.

Policy Gradient with MARL: To develop an algorithm for MARL, we would apply
the policy gradient theorem as mentioned in Theorem 3.1 in (Zhang et al.,
2018). For any $\theta\in\Theta$, let $\pi_{\theta}$ be a policy and
$J(\theta)$ denote the globally averaged return, we define the local advantage
function $A_{\theta}^{i}:\mathcal{S}\times\mathcal{A}$ $\rightarrow\mathbb{R}$
as
$A_{\theta}^{i}(s,a)=\mathcal{Q}_{\theta}(s,a)-\tilde{V}_{\theta}^{i}(s,a^{-1})$,
where $\mathcal{Q}_{\theta}$ and $A_{\theta}$ are the corresponding global
action-value and advantage functions and $\tilde{V}^{i}_{\theta}(s,a^{-1})$ is
the local value function. The gradient of the policy is given by

(3) $\nabla_{\theta^{i}}J(\theta)=\mathbb{E}_{s\sim
d_{\theta},a\sim\pi_{\theta}}[\nabla_{\theta^{i}}log\pi_{\theta^{i}}^{i}(s,a^{i})\cdot
A_{\theta}^{i}(s,a)]$

##  3\. Adversary Aware MARL Algorithm

###  3.1. Non Adversarial Setting

We first present the existing work under a non-adversarial setting. Here, he
actor-critic based consensus algorithm for the MARL setup is as follows: we
consider an agent specific local version of the global action-value function
$Q_{\theta}$ which we denote as $\mathcal{Q}(\omega^{i})$ where we hide the
state, action factors in the function which is implicit and
$\omega\in\mathbb{R}^{U}$, where the dimension $U$ is significantly less than
the joint state action space. In order to use the policy gradient theorem
discussed in the previous section to improve an agents’ policy, each agent
shares its local parameters $\omega^{i}$ with its neighbors on the network in
order to reach an estimate of $\mathcal{Q}_{\theta}$ that is consensual among
all agents in the network. Such distributed consensus algorithms have been
proposed earlier (Nedic and Ozdaglar, 2009) that guarantee convergence of the
local agent functions.

Algorithm 1 Adversary aware Networked actor-critic algorithm

1:Input: Initial values of the parameters $\mu 0^{i}$, $\omega_{0}^{i}$,
$\tilde{\omega}_{0}^{i}$, $\theta_{0}^{i}$, $\forall i\in\mathcal{N}$; the
initial state $s_{0}$ of the MDP and stepsizes $\\{\beta_{\omega,t}\\}_{t\geq
0}$ and $\\{\beta_{\theta,t}\\}_{t\geq 0},F,V(N)$

2:Repeat:

3: for all $i\in\mathcal{N}$ do

4: Observe state $s_{t+1}$ and reward $r^{i}_{t+1}$.

5: Update $\mu_{i+1}^{i}\leftarrow(1-\beta_{\omega,t})$ $\cdot$ $\mu_{t}^{i}$
\+ $\beta_{\omega,t}$ $\cdot$ $r_{t+1}^{i}$.

6: Select and execute action $a_{t+1}^{i}$ $\sim$
$\pi_{\theta_{t}^{i}}^{i}(s_{t+1},\cdot)$.

7: end for

8: Observe joint actions $a_{t+1}$ = $(a_{t+1}^{1},\ldots,a_{t+1}^{N})$.

9: for all $i\in\mathcal{N}$ do

10: Update $\delta_{t}^{i}\leftarrow r_{t+1}^{i}$ \- $\mu_{t}^{i}$ \+
$\mathcal{Q}_{t+1}(\omega_{t}^{i})$ \- $\mathcal{Q}_{t}(\omega_{t}^{i})$.

11: Critic Step: $\tilde{\omega}_{t}^{i}$ $\leftarrow$ $\omega_{t}^{i}$ \+
$\beta_{\omega,t}$ $\cdot$ $\delta_{t}^{i}$ $\cdot$
$\nabla_{\omega}\mathcal{Q}_{t}(\omega_{t}^{i})$.

12: Update $\mathcal{A}_{t}^{i}$ $\leftarrow$
$\mathcal{Q}_{t}(\omega_{t}^{i})$ \- $\sum_{a^{i}}\in\mathcal{A}^{i}$
$\pi_{\theta_{t}^{i}}(s_{t},a^{i}$ $\cdot$ $\mathcal{Q}(s_{t},a^{i},\\\
a^{-1};\omega_{t}^{i})$ .

13: Update $\psi_{t}^{i}$ $\leftarrow$ $\nabla_{\theta}^{i}$ log
$\pi_{\theta_{t}^{i}}^{i}(s_{t},a_{t}^{i})$.

14: Actor Step: $\theta_{t+1}^{i}$ $\leftarrow$ $\theta_{t}^{i}$ \+
$\beta_{\theta,t}$ $\cdot$ $\mathcal{A}_{t}^{i}$ $\cdot$ $\psi_{t}^{i}$.

15: Send $\omega_{t}^{i}$ to the neighbors
$\\{j\in\mathcal{N}:(i,j)\in\mathcal{E}_{t}\\}$ over the communication network
$\mathcal{G}_{t}$ .

16: end for

17: for all $i\in\mathcal{N}$ do

18: Gather neighbors $K_{i}$ from $V(N)$

19: Gather $\tau_{i}(t)$ by removing the highest and lowest $F$ states among
$K_{i}$

20: Consensus Step: $\omega_{t+1}^{i}$ $\leftarrow$
$\sum_{j\in\tau_{i}}c_{t}(i,j)$ $\cdot$ $\tilde{\omega}_{t}^{j}$.

21: end for

22: Update the iteration counter $t$ $\leftarrow$ $t+1$.

23:Until Convergence

###  3.2. Adversary Aware Decentralized MARL

In order to introduce the setup behind the MARL environment where malicious
agents are a part of the network, we first define certain notations based on
the setup of consensus based adversarial attacks in (Sundaram and Gharesifard,
2015). We consider an undirected network $\mathcal{N}={V(N),E(N)}$, where
$V(N)=\\{v1_{,}1\ldots,v_{n}\\}$ and $E(N)$ denotes the edges connecting pairs
of nodes. Note that $v_{i}$ is just an exaggerated notation to denote an agent
in a graph as opposed to $i$ that we have used in the previous sections but
they refer to the same agent. Denoting $\mathcal{K}_{i}$ as the neighborhood
vertices of $i\in V(N)$, for any $r\in\mathbb{N}$, a subset $S\subset V(N)$ is
said to be r-local, if $|\mathcal{K}_{i}\cap S|\leq r$ $\forall v_{i}\in
V(N)\setminus S$. That is to say, for r-local subset, there are at most $r$
nodes in the neighborhood of any vertex from $V(N)\setminus S$. For the setup
of adversarial attacks in the multi-agent system, we consider the following:
we consider that there are randomly chosen adversarial nodes in the network
such that for each node $v_{i}$, there cannot be more than $g$ fraction of its
neighbors who are adversaries, where $g\in[0,1)$. This fraction $g$ is known
to all nodes in the network, however the regular nodes do not know which or if
any of their neighbor nodes are adversaries. Here we assume that the
adversaries are restricted to form an F-local set , where $F$ is a non-
negative integer.

###  3.3. Algorithm

The overall actor-critic algorithm is detailed in Algorithm 1. In the actor
step, each agent, each agent improves it policy as shown in Line 14, where
$\beta_{\theta,t}>0$ is the stepsize. Note that both the actor and the critic
steps can be executed in a decentralized fashion without any centralized
training. For the consensus step, one important thing to note that since the
agents aim to optimize the globally averaged reward function $\overline{r}$,
the agents share their local parameters $\omega^{i}$ with their neighbors and
this allows the agents to improve their current policy using the policy
gradient theorem. The parameter sharing involves a consensus update using the
weight matrix $C_{t}=[c_{t}(i,j)]_{|\mathcal{N}|\times|\mathcal{N}|}$ such
that $c_{t}(i,j)$ is the weight on the message sent by agent $j$ to agent $i$
at time $t$. An important restriction in our model is that we only consider
$c_{t}(i,j)>0$ if agent $i$ and $j$ are neighbors of each other in the
network. Some discussions on the choice of $c_{t}(i,j)$ is in the later part
of this section.

We modify Algorithm 1 in (Zhang et al., 2018) to incorporate co-ordinated
responses by regular nodes in the presence of adversaries. At each time step,
the regular nodes $\\{v_{i}\\}$ gather $\tilde{\omega}_{t}^{j}$ for all
$v_{j}\in\mathcal{K}_{i}$ and remove the highest and lowest $F$ states and the
remaining nodes are denoted by $\tau_{i}(t)\in\mathcal{K}_{i}$. The consensus
action in the actor is then to aggregate the parameters from the neighbors in
$\tau_{i}(t)$. Note that, these updates are only done by the regular nodes and
the adversarial nodes aggregate the weights from the neighbors in any way they
wish. The dynamics of updates by the regular nodes are local and decentralized
since they do not require regular nodes to know anything beyond the signals
sent from their neighbors. As mentioned in (Sundaram and Gharesifard, 2015),
when the network $\mathcal{N}$ is time-invariant, the effective neighbor set
$\tau_{i}(t)$ is only a function of the states of neighbors of $v_{i}$ at time
step $t$. This filtering is also closely connected to bandit based ranking of
arms where the distribution of the means of the arms decide which arm would be
picked, albeit here instead of picking one arm, the agent selects one or
multiple arms (Even-Dar et al., 2006). The rest of the consensus style
algorithm in this decentralized environment is the same as Algorithm 1 in
(Zhang et al., 2018).

The standard assumptions for the consensus matrix is defined in Assumption 4.4
of (Zhang et al., 2018) and while there can be many ways to define the weight,
one popular way is to consider the notion of Metropolis weights (Xiao et al.,
2005)

$\displaystyle c_{t}(i,j)$
$\displaystyle=\Big{\\{}1+max[k_{t}(i),k_{t}(j)]\Big{\\}}^{-1},\ \
\forall(i,j)\in\mathcal{E}_{t}$ $\displaystyle c_{t}(i,j)$
$\displaystyle=1-\sum_{j\in\tau_{i}}c_{t}(i,j)\forall i\in\mathcal{N}$

where $\mathcal{E}_{t}$ denotes the time varying edges and $k_{t}(i)$ is the
degree of the agent $i$ in the time varying network.

##  4\. Related Work

One of the assumptions in the process of social decision making is that
individuals following a learning trajectory (despite each individual having a
limited memory) successfully converges to the best option for the collective
population. While individuals participate in decision making where they
associate different risk and rewards in an uncertain environment, they also
tend to incorporate beliefs from their immediate neighbors in a networked
environment, a phenomenon that has played a critical role in co-operative
multi-agent systems (Celis et al., 2017). This begets the question as to how
resilient are these consensus algorithms and how can agents adapt to
decentralized training in the presence of malicious agents or adversaries. We
briefly highly some notable studies done previously at the intersection of
MARL and

MARL: The field of MARL has evolved very rapidly over the past few years. The
collective goal of the multi-agent system is to either reach a stable and
consensus state for all agents (Fax and Murray, 2004), or solve a static
optimization problem in a distributed fashion(Nedic and Ozdaglar, 2009).
However, competition and collaboration always emerge between autonomous agents
that learn by reinforcement over finite horizons. Some of the common
approaches on modeling and solving cooperative multi-agent reinforcement
learning problems include: (1) independent learners (Foerster et al., 2017;
Laurent et al., 2011), (2) fully observable critic (Mao et al., 2018), (3)
value function factorization (Son et al., 2019), (4) consensus based RL
(Cassano et al., 2020), and (5) learning to communicate (Mordatch and Abbeel,
2018).

Decentralized Networked MARL: The challenge with fully co-operative multi-
agent systems is that fully cooperative systems (Dec-POMDPs) are significantly
harder to solve than single agent RL problems due to the combinatorial
explosion in the joint action state spaces combining all agents. For this
reason CTDE have become popular. The prevailing MARL paradigm of centralised
training with decentralised execution (CTDE) [25,29,21] assumes a training
stage during which the learning algorithm can access data from all agents to
learn decentralised (locally-executable) agent policies. CTDE algorithms such
as COMA (Foerster et al., 2018) learn powerful critic networks conditioned on
joint observations and actions of all agents. Other extensions of MADDP
include shared experience actor-critic frameworks (Christianos et al., 2020)
for efficient MARL exploration by sharing experience amongst agents as opposed
to MADDPG which only reinforces an agent’s own explored actions. There have
also been studies conducted on making these RL systems in networked systems
scalabale (Qu et al., 2020).

Decision making for social systems: There are two broad avenues of research on
decision making in interactive environments considering social systems. The
first group of studies focus on game theoretic environments where agents have
similar adaptation and learning abilities and so the actions of each agent
affect the task achievement of the other agents (Bowling and Veloso, 2000;
Nakayama et al., 2017; Reverdy et al., 2014). The payoff of the agents in
these games depend on whether they are purely collaborative or competitive or
a mix. In a recent study on combining Reinforcement Learning with Agent Based
Modeling (Sert et al., 2020), the authors address the self‑organizing dynamics
of social segregation and explore the space of possibilities that emerge from
considering different types of rewards. The second group of studies focuses on
social influence from other agents as the intrinsic factor for the decision
making. Social influence has been known to be an intrinsic factor in agents’
choices (Sarkar et al., 2020; Schlag, 1998), and in recent studies, there have
been attempts to propose a unified mechanism for coordination in MARL by
rewarding agents for having causal influence over other agents’ actions
(Jaques et al., 2019; Xie et al., 2021). Adversarial attacks during training
can potentially impact the choices made by agents and impact the consensus
algorithms as studied in (Figura et al., 2021) and this is important to solve
since attacks through social influence are a common tool for adversaries
(Sarkar et al., 2019).

##  5\. Conclusion

In this work, we discussed a simple strategy for agents to mitigate the
nuances of adversarial agents in fully decentralized MARL where the agents are
connected via a time varying network. We relied on the algorithm of (Zhang et
al., 2018) and modified the actor-critic algorithm in the presence of one or
more adversaries. We did not provide a formal argument for the convergence of
the consensus updates with our adversarial framework which we leave as one of
the immediate next steps of this work. A second direction where this kind of
work can be adapted is to use attention based actor critic algorithms
(iqbal2019actor) that modify the action value function by parameterizing it
with agent specific attention weights. We can then use reward values to
optimize for the attention weights over the time varying network.

## References

  * (1)
  * Bowling and Veloso (2000) Michael Bowling and Manuela Veloso. 2000.  _An analysis of stochastic game theory for multiagent reinforcement learning_.  Technical Report. Carnegie-Mellon Univ Pittsburgh Pa School of Computer Science.
  * Busoniu et al. (2008) Lucian Busoniu, Robert Babuska, and Bart De Schutter. 2008\.  A comprehensive survey of multiagent reinforcement learning.  _IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)_ 38, 2 (2008), 156–172.
  * Cassano et al. (2020) Lucas Cassano, Kun Yuan, and Ali H Sayed. 2020.  Multiagent fully decentralized value function learning with linear convergence rates.  _IEEE Trans. Automat. Control_ 66, 4 (2020), 1497–1512.
  * Celis et al. (2017) L Elisa Celis, Peter M Krafft, and Nisheeth K Vishnoi. 2017\.  A distributed learning dynamics in social groups. In _Proceedings of the ACM Symposium on Principles of Distributed Computing_. 441–450.
  * Christianos et al. (2020) Filippos Christianos, Lukas Schäfer, and Stefano Albrecht. 2020.  Shared experience actor-critic for multi-agent reinforcement learning.  _Advances in neural information processing systems_ 33 (2020), 10707–10717.
  * Even-Dar et al. (2006) Eyal Even-Dar, Shie Mannor, Yishay Mansour, and Sridhar Mahadevan. 2006.  Action Elimination and Stopping Conditions for the Multi-Armed Bandit and Reinforcement Learning Problems.  _Journal of machine learning research_ 7, 6 (2006).
  * Fax and Murray (2004) J Alexander Fax and Richard M Murray. 2004.  Information flow and cooperative control of vehicle formations.  _IEEE transactions on automatic control_ 49, 9 (2004), 1465–1476.
  * Figura et al. (2021) Martin Figura, Yixuan Lin, Ji Liu, and Vijay Gupta. 2021\.  Resilient Consensus-based Multi-agent Reinforcement Learning.  _arXiv preprint arXiv:2111.06776_ (2021).
  * Foerster et al. (2018) Jakob Foerster, Gregory Farquhar, Triantafyllos Afouras, Nantas Nardelli, and Shimon Whiteson. 2018\.  Counterfactual multi-agent policy gradients. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 32.
  * Foerster et al. (2017) Jakob Foerster, Nantas Nardelli, Gregory Farquhar, Triantafyllos Afouras, Philip HS Torr, Pushmeet Kohli, and Shimon Whiteson. 2017.  Stabilising experience replay for deep multi-agent reinforcement learning. In _International conference on machine learning_. PMLR, 1146–1155.
  * Jaques et al. (2019) Natasha Jaques, Angeliki Lazaridou, Edward Hughes, Caglar Gulcehre, Pedro Ortega, DJ Strouse, Joel Z Leibo, and Nando De Freitas. 2019\.  Social influence as intrinsic motivation for multi-agent deep reinforcement learning. In _International conference on machine learning_. PMLR, 3040–3049.
  * Kok and Vlassis (2006) Jelle R Kok and Nikos Vlassis. 2006.  Collaborative multiagent reinforcement learning by payoff propagation.  _Journal of Machine Learning Research_ 7 (2006), 1789–1828.
  * Laurent et al. (2011) Guillaume J Laurent, Laëtitia Matignon, Le Fort-Piat, et al. 2011\.  The world of independent learners is not Markovian.  _International Journal of Knowledge-based and Intelligent Engineering Systems_ 15, 1 (2011), 55–64.
  * Lowe et al. (2017) Ryan Lowe, Yi I Wu, Aviv Tamar, Jean Harb, OpenAI Pieter Abbeel, and Igor Mordatch. 2017.  Multi-agent actor-critic for mixed cooperative-competitive environments.  _Advances in neural information processing systems_ 30 (2017).
  * Lyu et al. (2021) Xueguang Lyu, Yuchen Xiao, Brett Daley, and Christopher Amato. 2021.  Contrasting centralized and decentralized critics in multi-agent reinforcement learning.  _arXiv preprint arXiv:2102.04402_ (2021).
  * Mao et al. (2018) Hangyu Mao, Zhengchao Zhang, Zhen Xiao, and Zhibo Gong. 2018\.  Modelling the dynamic joint policy of teammates with attention multi-agent DDPG.  _arXiv preprint arXiv:1811.07029_ (2018).
  * Mordatch and Abbeel (2018) Igor Mordatch and Pieter Abbeel. 2018.  Emergence of grounded compositional language in multi-agent populations. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 32.
  * Nakayama et al. (2017) Kazuaki Nakayama, Masato Hisakado, and Shintaro Mori. 2017\.  Nash equilibrium of social-learning agents in a restless multiarmed bandit game.  _Scientific reports_ 7, 1 (2017), 1937.
  * Nedic and Ozdaglar (2009) Angelia Nedic and Asuman Ozdaglar. 2009.  Distributed subgradient methods for multi-agent optimization.  _IEEE Trans. Automat. Control_ 54, 1 (2009), 48–61.
  * Qu et al. (2020) Guannan Qu, Yiheng Lin, Adam Wierman, and Na Li. 2020\.  Scalable multi-agent reinforcement learning for networked systems with average reward.  _Advances in Neural Information Processing Systems_ 33 (2020), 2074–2086.
  * Reverdy et al. (2014) Paul B Reverdy, Vaibhav Srivastava, and Naomi Ehrich Leonard. 2014\.  Modeling human decision making in generalized Gaussian multiarmed bandits.  _Proc. IEEE_ 102, 4 (2014), 544–571.
  * Sarkar et al. (2019) Soumajyoti Sarkar, Paulo Shakarian, Mika Armenta, Danielle Sanchez, and Kiran Lakkaraju. 2019\.  Can social influence be exploited to compromise security: An online experimental evaluation. In _Proceedings of the 2019 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining_. 593–596.
  * Sarkar et al. (2020) Soumajyoti Sarkar, Paulo Shakarian, Danielle Sanchez, Mika Armenta, and Kiran Lakkaraju. 2020\.  Use of a controlled experiment and computational models to measure the impact of sequential peer exposures on decision making.  _Plos one_ 15, 7 (2020), e0234875.
  * Schlag (1998) Karl H Schlag. 1998\.  Why imitate, and if so, how?: A boundedly rational approach to multi-armed bandits.  _Journal of economic theory_ 78, 1 (1998), 130–156.
  * Sert et al. (2020) Egemen Sert, Yaneer Bar-Yam, and Alfredo J Morales. 2020\.  Segregation dynamics with reinforcement learning and agent based modeling.  _Scientific reports_ 10, 1 (2020), 11771.
  * Son et al. (2019) Kyunghwan Son, Daewoo Kim, Wan Ju Kang, David Earl Hostallero, and Yung Yi. 2019.  Qtran: Learning to factorize with transformation for cooperative multi-agent reinforcement learning. In _International conference on machine learning_. PMLR, 5887–5896.
  * Sundaram and Gharesifard (2015) Shreyas Sundaram and Bahman Gharesifard. 2015.  Consensus-based distributed optimization with malicious nodes. In _2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton)_. IEEE, 244–249.
  * Willemsen et al. (2021) Daniël Willemsen, Mario Coppola, and Guido CHE de Croon. 2021\.  MAMBPO: Sample-efficient multi-robot reinforcement learning using learned world models. In _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 5635–5640.
  * Xiao et al. (2005) Lin Xiao, Stephen Boyd, and Sanjay Lall. 2005.  A scheme for robust distributed sensor fusion based on average consensus. In _IPSN 2005. Fourth International Symposium on Information Processing in Sensor Networks, 2005._ IEEE, 63–70.
  * Xie et al. (2021) Annie Xie, Dylan Losey, Ryan Tolsma, Chelsea Finn, and Dorsa Sadigh. 2021.  Learning latent representations to influence multi-agent interaction. In _Conference on robot learning_. PMLR, 575–588.
  * Zhang et al. (2018) Kaiqing Zhang, Zhuoran Yang, Han Liu, Tong Zhang, and Tamer Basar. 2018.  Fully decentralized multi-agent reinforcement learning with networked agents. In _International Conference on Machine Learning_. PMLR, 5872–5881.


Pranshu Gupta

Indian Institute of Technology, Kanpur

Algorithms Inspired by Nature: A Survey

###### Contents

  1. 1 Introduction
  2. 2 Ant Colony Optimization
    1. 2.1 Ant Colony Optimization Metaheuristic
    2. 2.2 Applications of Ant Algorithms
      1. 2.2.1 Job Shop Scheduling
      2. 2.2.2 Aircraft Conflict Resolution
      3. 2.2.3 Ant Based Clustering
  3. 3 Simulated Annealing
    1. 3.1 Solving Job Shop Scheduling with Simulated Annealing
  4. 4 Particle Swarm Optimization
    1. 4.1 Applications of Particle Swarm Optimization
      1. 4.1.1 Solving Job Shop Scheduling with Particle Swarm Optimization
      2. 4.1.2 A Hybrid of SA and PSO to solve JSSP
  5. 5 Evolutionary Algorithms
    1. 5.1 Encoding a solution as a Chromosome
      1. 5.1.1 Binary Coding
      2. 5.1.2 Delta Coding
      3. 5.1.3 Fuzzy Coding
    2. 5.2 Selection and Replacement
      1. 5.2.1 Roulette Wheel Selection
      2. 5.2.2 Ranking Selection
      3. 5.2.3 Tournament Selection
      4. 5.2.4 Replacement
    3. 5.3 Reproduction
    4. 5.4 Applications of Evolutionary Algorithms
      1. 5.4.1 Symbolic Regression with Genetic Algorithms
      2. 5.4.2 Support Vector Regression with Genetic Algorithms
  6. 6 Conclusion

##  Introduction

Nature is known to be the best optimizer. Natural processes most often than
not reach an optimal equilibrium. Scientists have always strived to understand
and model such processes. Thus, many algorithms exist today that are inspired
by nature. Many of these algorithms and heuristics can be used to solve
problems for which no polynomial time algorithms exist, such as Job Shop
Scheduling and many other Combinatorial Optimization problems. We will discuss
some of these algorithms and heuristics and how they help us solve complex
problems of practical importance. We will focus mainly on Ant Colony
Optimization, Particle Swarm Optimization, Simulated Annealing and
Evolutionary Algorithms.

######  Remark 0.1.

A heuristic can be thought of as a rule of thumb that will hopefully find a
good solution for a given problem. These techniques are mainly applied to
complex convex multivariate constrained combinatorial optimization problems
which are usually NP complete or harder.

######  Remark 0.2.

A metaheuristic is a higher-level procedure or heuristic designed to find,
generate, or select a heuristic that may provide a sufficiently good solution
to an optimization problem, especially with incomplete or imperfect
information or limited computation capacity. [Bianchi et al. (2009)]

Ant Colony Optimization utilizes a family of heuristics inspired by the
foraging behavior of Ants. Ants are highly co-operative social insects, they
search for food in huge groups and communicate indirectly using a chemical
called pheromone. These heuristics model individual ants as independent
solution finders roaming in the solution space of a problem and communicating
via a shared variable(s) i.e. the pheromone.

When a liquid is cooled slowly, its particles arrange themselves in some
configuration that allows minimum energy state - this process is called
annealing. Simulated Annealing tries to model this phenomenon by drawing an
analogy between between the behavior of multi-body systems (e.g. liquids) in
thermal equilibrium at a finite temperature and the behaviour of the objective
function of a combinatorial optimization problem at different values of the
variables involved. [Kirkpatrick et al. (1983)]

Particle Swarm Optimization is a another optimization heuristic. It relies on
a swarm of particles that are initialized with some configuration in the
solution space of a given problem. These particles then explore the solution
space individually while communicating with other particles around them to
find the best solution. After some iterations these particles are expected to
make swarms around optimal solution(s). PSO relies on randomization in the
movement of particles to avoid getting stuck in local optima.

Evolutionary Algorithms are a family of algorithms inspired from Neo-Darwinism
paradigm of biological evolution.

$\text{Neo-Darwininsm}=\text{Darwinism}+\text{Selectionism}+\text{Mendelian
Genetics}$

These algorithms model solutions as individual species with their genetic code
representing the solution itself. These individuals later reproduce to create
new individuals with different genetic code. Reproduction involves choosing
healthy parents and then crossing over their genes to create an offspring. The
optimality of the solution instance represented by the genetic code of the
offspring determines its fitness. Then the notion of ”survival of the fittest”
comes into play. The best individuals survive and then they reproduce for the
next generational cycle. Thus, this model utilizes both exploration and
exploitation of the solution space.

Genetic Algorithms, Evolutionary Strategies, Differential Evolution etc., are
a few variants of EAs.

##  Ant Colony Optimization

Ant Algorithms are a family of algorithms inspired by the behavior of real ant
colonies. Ants are social insects and they always work together to ensure the
survival of the colony as a whole. When in search of food, ants deposit a
chemical substance called pheromone on their path. If they succeed in finding
a food source they again deposit pheromone on their path back to the colony.
The amount of pheromone they deposit on the path often is determined by the
quality and quantity of the food source that the path leads to. Also, if the
multiple ants find the same food source, the one which followed the shorter
path will deposit the pheromone first. When looking for food the ants also
consider the pheromone deposited by other ants to aid their own search. Thus,
the shorter paths and the paths which lead to better food sources will be
followed by more and more ants and the pheromone on them will be further
enriched (positive feedback in search). Also, because pheromone evaporates in
atmosphere with the passage of time, paths which are no longer useful or
followed will lose the pheromone deposited on them. However, ants are always
free to choose their own paths ensuring the explore for new food sources as
well instead of just clinging upon only one.

This behavior of ants can be utilized to design algorithms to solve search and
optimization problems. Dorigo et al. (1999) were the first to model foraging
behavior of ant colonies and design a metaheuristic that can be used to solve
problems of practical importance.

###  Ant Colony Optimization Metaheuristic

The ACO metaheuristic is usually defined in the spirit of Combinatorial
Optimization problems. It is indeed one of the family of approximate
algorithms that are used to solve hard CO problems in reasonable amount of
time. Let’s define a CO problem.

######  Definition 0.1.

A Combinatorial Optimization problem $P$ is defined as $P=(S,\omega,f)$ where
$S$ is the solution space for the given problem, $\omega$ is the constraint
set on the solution parameters and $f:S\rightarrow{\mathbb{R}}^{+}$ is the
objective function to be minimized.

The solution set can be defined as a set of feasible solutions. A feasible
solution is a set of $N$ discrete variables $X_{i}\in D_{i}$ where
$i\in{1,...n}$ which satisfy the solution constraints given by the problem
instance.

The pheromone trail model is the is the most crucial component of any ACO
framework. It is used to probabilistically build the solution components of
the problem instance. A solution component $c_{ij}$ is defined as the
assignment of a variable $X_{i}$ to a value $v_{j}\in D_{i}$. The pheromone
trail corresponding to this solution component is denoted by $T_{ij}$ with its
value being $\tau_{ij}$. The set of all solution components is denoted by $C$
and the vector of pheromone trails is denoted by $T$.

For ACO, we define a complete graph which has a node corresponding to each
possible solution component called the Construction Graph. In each iteration
of ACO, each ant moves from one component to another effectively assigning
values to the discrete variables (taking care of the constraints and
considering deposited pheromone) and thus creating a feasible solution. At the
end of every iteration, each ant would have constructed a feasible solution.
Then, we compute the fitness of these solutions and deposit pheromone trail
values on the edges that the ant which created the solution followed.

1 Initialize Pheromone Trail values $T$

2 The best solution $S_{bs}\leftarrow$ NULL

3 while  _termination conditions not met_ do

4 Set of solutions found in this iteration $S_{iter}\leftarrow\phi$

5 foreach  _ant in colony_ do

6 $s\leftarrow$ ConstructSolution($T$)

7 if  _$s$ is feasible_ then

8 $s\leftarrow$ LocalSearch($s$)

9 if  _$f(s) <f(s_{bs})$ OR $s_{bs}$ is NULL_ then

10 $s_{bs}\leftarrow s$

11

12 end if

13 $S_{iter}\leftarrow S_{iter}\cup\\{s\\}$

14

15 end if

16

17 end foreach

18 ApplyPheromoneUpdate($T,S_{iter},s_{bs}$)

19

20 end while

21Return $s_{bs}$

Algorithm 1 Generic ACO Algorithm [Dorigo and Blum (2005)]

The probability that an ant will choose a solution component $c_{ij}$ as part
of its solution conditioned on the set of solution components it has already
chosen (the partial solution $s^{p}$) is given by:

$p(c_{ij}|s^{p})=\frac{\tau_{ij}^{\alpha}\cdot c_{ij}^{\beta}}{\sum_{c_{lk}\in
N(s^{p})}\tau_{lk}^{\alpha}\cdot c_{lk}^{\beta}}$

Here, $c_{ij}$ denotes the problem specific cost of choosing the assignment
$j$ for decision variable $i$. $\alpha$ and $\beta$ are parameters than are
used to tune the relative importance of the problem specific costs and
pheromone trail value corresponding to solution component. $N(s^{p})$ is the
set of allowed solution components given that the partial solution $s^{p}$ has
already been constructed. $N(s^{p})$ depends on the constraint set of the
problem instance. If the constraint set does not allow the selection of a
solution component given the partial solution, the corresponding transition
probability is 0.

The pheromone update rule is given as follows

$\tau_{ij}\leftarrow(1-\rho)\tau_{ij}+Q\sum_{k=1}^{m}\frac{1}{f(s_{k})}$

Here, $\rho$ is the evaporation coefficient which models the atmospheric
evaporation of the pheromone, Q is some constant, $s_{k}$ is the solution
constructed by the $k^{th}$ ant.

There are other variants of ACO which use different pheromone update rules.
Consider the MIN-MAX ACO, the pheromone update rule for it is given by the
following equation where $s_{iterbs}$ is the best solution found in this
iteration. [Dorigo et al. (2006)]

$\tau_{ij}\leftarrow(1-\rho)\tau_{ij}+\frac{1}{f(s_{iterbs})}$

In the Ant Colony System Framework, instantaneous updates are also made to the
pheromone trails. This slightly decreases the pheromone value allowing other
ants to possible find diverse solutions instead of premature convergence to a
sub-optimal solution.

$\tau_{ij}\leftarrow(1-\phi)\tau_{ij}+\phi\tau_{0}$

###  Applications of Ant Algorithms

####  Job Shop Scheduling

A Job Shop Scheduling problem consists of a set of $n$ Jobs and a set of $m$
machines. Each job consists of a sequence operations each to be processed on a
specific machine. Now, because the resources have to be shared and each
operation takes a given amount of time, we want to minimize the total time
span required to finish all the jobs - the makespan. There are many techniques
to solve this problem - Shifting Bottleneck, SA, ACO, PSO a well as EAs can be
used to solve JSSP.

Every JSSP instance can be represented by a disjunctive graph [See Figure 1]
in which we have one node for each operation with directed edges from $o_{i}$
to $o_{j}$ if the latter comes after the former in the constraint sequence
given in the problem. These constraint sequences are generally given for
operations that belong to the same job. We dont have a complete ordering
across all the operations from different jobs. So, in the disjunctive graph we
have un-directed edges between the operations belonging to different jobs but
the same machine. If we define a machine specific ordering of operations on
each machine, we have essentially solved the problem. We also have a start and
finish node to denote the start and finish of the complete sequence of
operations [Yamada and Nakano (1997)]

Figure 1: Disjunctive Graph for a JSSP with 3 machines and 3 Jobs each with 3
operations. $O_{ij}$ denotes an operation belonging to job $i$ which has to
processed on machine $j$, $p_{ij}$ denotes the processing time corresponding
to the operation. Image taken from [Yamada and Nakano (1997)]

In ACO, we create the construction graph for ant traversal using the
disjunctive graph of the problem instance. We add two directed edges in the
construction graph for each undirected edge in the disjunctive graph. An ant
will choose on of those edges and hence the machine specific ordering will be
found. Now, an ant walk from start node to the finish node in the construction
graph will essentially give a complete solution. Note that the ant visits all
the nodes in the graph one by one. The pheromone update rule captures the
makespan as the fitness function and the transition probabilities respect the
feasibility of the solution by making the probability zero if the transition
creates a cycle in the ant walk, as a cyclic order is non-realizable. [Blum
and Sampels (2002)]

####  Aircraft Conflict Resolution

The Aircraft Conflict Resolution problem is of great practical importance. It
is combinatorial optimization problem in which a set of $n$ aircrafts
approaching a single point have to deviate from their original course in order
to avoid collision.

Figure 2: Aircraft Conflict Problem. Image taken from [Durand and Alliot
(2009)]

Durand and Alliot (2009) have solved a simplified version of this problem in
which all the aircrafts are approaching with the same speed, time is
discretized, the set of possible manoeuvres is finte (say
$10^{o},20^{o},30^{o}$) and only two manoeuvres are to be applied to each
aircraft. Then this becomes a combinatorial optimization problem wherein we
have to choose the manoeuvres at 2 time steps. The construction graph they
have designed is shown in Figure 3. Each node is present at a discrete time
step and there are three kinds of nodes U, V and W. The only allowed
transition between differently typed nodes are U to V and V to W.

Figure 3: Construction Graph for Aircraft Conflict Problem. Image taken from
[Durand and Alliot (2009)]

Now, in each iteration of ACO, a set of $n$ ants create a single solution for
the problem. Each of these ants corresponds to an aircraft and the walk done
by the the ant defines the manoeuvres for the aircraft. We take this set on
$n$ walks and see if they give a conflict free plan, if they do, the pheromone
trail is updated according to the delay due to the course correction as the
fitness function. They have demonstrated that this technique can resolve
conflicts for upto 30 aircrafts in reasonable time.

####  Ant Based Clustering

Some species of ants are known to arrange their food supplies in clusters on
the basis of different properties of the food items (size etc.). We can also
model this behavior of ants to solve clustering problems as well.

Handl et al. (2006) have proposed a simple ant based clustering algorithm that
gives natural clusters. Unlike K-Means and other traditional clustering
techniques which often need the programmer to specify the number of cluster
centers beforehand, this technique gives natural clusters all by itself. They
have experimentally shown that it is very robust and gives clusters that are
inherent from the dataset, even if the clusters overlap.

1 Randomly scatter all the data points in a toroidal grid

2 foreach  _ant in colony_ do

3 Ant picks a randomly chosen datapoint ($i$)

4 Ant places itself at a randomly selected location in the grid

5

6 end foreach

7while  _termination conditions not met_ do

8 Choose a random ant and displace it by a given stepsize

9 Drop the datapoint ($i$) at this location in the grid with probability given
by  $p=\bigg{(}\frac{f^{*}(i)}{0.3+f(i)}\bigg{)}^{2}$ if  _datapoint was
dropped_ then

10 while  _another datapoint is not picked_ do

11 Randomly choose a datapoint ($i$)

12 Pick it up with the probability given by
$p=\bigg{(}\frac{0.1}{0.3+f(i)}\bigg{)}^{2}$

13 end while

14

15 end if

16

17 end while

Algorithm 2 Ant Based Clustering [Handl et al. (2006)]

The function $f(i)$ is the Lumer and Faieta neighborhood function defined as:

$f(i)=max\bigg{(}0,\frac{1}{\sigma^{2}}\sum_{j\in
L}\bigg{(}1-\frac{\delta(i,j)}{\alpha}\bigg{)}\bigg{)}$

where $\sigma^{2}$ is the size of local neighborhood in the toroidal grid,
$\delta(i,j)$ is the distance function for the datapoints to be clustered, $L$
is the set of datapoints currently inside the $\sigma^{2}$ neighborhood in the
grid. $\alpha$ is data specific scaling parameter. The modified LF function
$f^{*}$ is given as follows:

$f^{*}(i)=\begin{cases}f(i)&\quad\text{if $\forall j$
$\big{(}1-\frac{\delta(i,j)}{\alpha}\big{)}$}\\\
$0$&\quad\text{otherwise}\end{cases}$

This modified LF neighborhood function inflicts heavy penalties for high
dissimilarities thus improving the spatial separation between the clusters.

##  Simulated Annealing

Annealing is the process of slowly cooling molten solids to get them in a
crystallized state. If the object is cooled rapidly, its particle do not reach
an equilibrium state and the crystals formed are small, thus reducing the
overall strength of the object. Therefore, the cooling plan in this procedure
is very important. In optimization settings, the simulated annealing algorithm
we model the particles of the cooling object as the parameters of a possible
solution of the given problem. The temperature is modeled as a variable which
defines the transition probability of the selected solution. Each run of the
algorithm consists of multiple iterations until an equilibrium (convergence)
is reached. After that we change the transition probability of the solution by
lowering the temperature. We stop the runs when the lower bound of temperature
reached. See algorithm 3

1 Choose a random configuration $X_{i}$, select the initial system temperture
$T_{0}$ and decide the cooling schedule (usually exponential)

2 Evaluate Cost $C(X_{i})$

3 while  _temperature ¿ lower bound_ do

4 Perturb $X_{i}$ to get a neighbor configuration $X_{i+1}$

5 Evaluate Cost $C(X_{i+1})$

6 if  _$C(X_{i+1}) <C(X_{i})$_ then

7 Set the new solution configuration to $C(X_{i+1})$

8

9 else

10 Set the new solution configuration to $C(X_{i+1})$ with probability p
$p=e^{-(C(X_{i+1})-C(X_{i}))/T_{k}}$

11 end if

12 if  _Equilibrium reached at $T$_ then

13 Reduce the system temperature to  $T_{k+1}=(\frac{T_{1}}{T_{0}})^{k}T_{k}$

14 end if

15

16 end while

Algorithm 3 Simulated Annealing

###  Solving Job Shop Scheduling with Simulated Annealing

The objective of the Job Shop Scheduling problem is to find out the optimal
sequence of jobs on a given set of machines so that the makespan is minimized.
JSSP is a NP Hard problem and Simulated Annealing can be used to solve it. Van
Laarhoven et al. (1992) have shown that SA performs very well on JSSP
instances. They modeled the configuration as the set of machine specific Job
sequences with the makespan time as their costs. The neighborhood structure is
defined as the set of configurations that can be obtained by reversing the
order of two successive operations (on the same machine) on a critical path in
the disjunctive graph. The perturbations in the SA give one of the neighbor
solutions. This way the SA procedure can be adopted to solve any JSSP
instance. Similarly, SA can be used to other combinatorial optimization
problems like TSP and Physical Chip Design [See Kirkpatrick et al. (1983)]

##  Particle Swarm Optimization

In Particle Swarm Optimization, a number of particles are initially scattered
in the search space of a given problem. The fitness of each particle can be
evaluated using its location in the search space. Now, in each iteration every
particle displaces itself to a new position. This new position is determined
by the history of best positions of this particle as well as some other
particles in the swarm. It is in fact the weighted mean of the particle’s own
best position and the neighbors best position. The weights are picked
uniformly from range of positive numbers. Eventually, the swarm as a whole is
likely to move close to the optimum of the fitness function. See Algorithm 4
for a generic PSO technique.

1 Initialize a population of particles at random positions in the search space

2 while  _termination conditions not met_ do

3 For each particle $i$ evaluate the fitness function at its current position

4 Update the best fitness value $pbest_{i}$ and position $p_{i}$ for this
particle as per step 3

5 Identify the best fitness particle $g$ in the neighborhood

6 Update the position and velocity of the particle $i$ as per the following
equations  $\displaystyle v_{i}$ $\displaystyle\leftarrow
v_{i}+U(0,\phi_{1})\otimes(p_{i}-x_{i})+U(0,\phi_{2})\otimes(p_{g}-x_{i})$
$\displaystyle x_{i}$ $\displaystyle\leftarrow x_{i}+v_{i}$

7 end while

Algorithm 4 Particle Swarm Optimization

This works quite well if the solution space of the problem is continuous,
unlike the discrete valued solution spaces of problems like TSP, JSP etc.
Defining velocities is a crucial part of the procedure, but according to Clerc
(2010), where a version of PSO for TSP is presented, velocity can be viewed as
an operator which when applied to the current position gives a new position
from the neighborhood.

There are many other variants of the Standard PSO, one in which the velocity
of the particles is bounded, another in which the updates are determined not
only by the particle’s own history and its best neighbor but its entire
neighborhood - Fully Informed Particle Swarm Optimization [Poli et al.
(2007)]. FIPSO has been demonstrated to work better than the standard PSO -
but it is sensitive to the way the neighborhood is defined.

###  Applications of Particle Swarm Optimization

####  Solving Job Shop Scheduling with Particle Swarm Optimization

In PSO, the configuration and neighborhood definitions are the same as those
in Simulated Annealing. However, we need to define the notion of velocity, the
distance between two configurations and the sum of two velocities. The
velocity itself can be defined as a permutation function that gives a new
machine specific sequence of operations from the old one. Thus, the velocity
can be defined as a list of pair of operations that belong to same machine,
these pairs are to be exchanged in the machine specific ordering to get the
new position i.e. a list of transposition. The distance between two positions
is the velocity that must be applied to the first position to reach the
second. Finally, the sum of two velocities is the order preserving truncated
union of the transposition sets used to define the two velocities. With these
definitions at hand we can go ahead and use PSO to solve JSSP. See Shao et al.
(2013)

####  A Hybrid of SA and PSO to solve JSSP

In a hybrid solution proposed by Shao et al. (2013), we initialize a large
number of particles using an approach proposed by Pezzella et al. (2008). The
new position for each particle is found by simulated annealing approach. Thus,
the structure of the overall algorithm remains like that of PSO (with SA
embedded in it).

The advantage of this hybrid framework is that it takes the best from SA -
local search and PSO - global search. Thus, avoiding local optima better than
traditional SA, maintaining diversity in solutions and identifying better
neighbors than the traditional PSO.

##  Evolutionary Algorithms

Evolutionary Algorithms are a class of general purpose algorithms that are
designed to solve optimization problems that are large, complex,
discontinuous, non-differentiable, multi-modal and possibly constrained as
well. EAs are highly parallelizable, simple and are known to almost always
reach a near optimum or the global optimum.

Any Evolutionary Algorithm generally comprises of the following components - a
population generator, a population selector, fitness estimator and
reproduction operators (crossover and mutation). The pseudocode of a generic
EA can be given as Algorithm 5

1 $t=0$

2 Randomly generate initial population $P(0)$

3 while  _termination condition not met_ do

4 Evaluate fitness of each individual of $P(t)$

5 Select individuals as parents from $P(t)$ based on their fitness

6 Apply search operators (crossover and mutation) on parents and generate
$P(t+1)$

7 $t:=t+1$

8

9 end while

Algorithm 5 A Generic Evolutionary Algorithm [Du and Swamy (2016)]

$P(t)$ denotes the current set of individuals in the population. Each
individual in the population is called a chromosome which represents a
solution to the given problem. The chromosome is denoted as a string of
elements called genes which represent individual parameters of the solution.
An individual’s chromosome is termed as the genotype and genotype when denoted
along with fitness and other properties of the corresponding solution is
called a phenotype.

###  Encoding a solution as a Chromosome

There are different techniques to represent a solution as a chromosome, each
having their pros and cons. The performance of these algorithms depends on the
encoding technique used, because the crossover and mutation operators are
applied on these representations only to explore the solution space. Some of
the notable representation techniques are Binary Coding, Gray Coding, Delta
Coding and Fuzzy Coding.

####  Binary Coding

Genetic Algorithms use Binary Encoding wherein each chromosome is an array of
genes $x$. Each of these genes $x_{i}$ represents some parameter of the
solution in form of binary string of length $l$. To recover the actual value
of the parameter whose value lies in the interval $[x_{i}^{-},x_{i}^{+}]$ we
use Equation 1 [Du and Swamy (2016)]

$x_{i}=x_{i}^{-}+(x_{i}^{+}-x_{i}^{-})\frac{1}{2^{l_{i}}-1}(\sum_{j=0}^{l_{i}-1}s_{j}2^{i})$
(1)

Binary Coding is a fixed length and static coding which sometimes leads to
slower convergence and lower accuracy when it reaches near the optimal
solution. Therefore, new dynamic coding techniques were developed, such as
Delta and Fuzzy Coding.

####  Delta Coding

Delta Coding is a dynamic coding technique in which the representation of the
search space changes in after each iteration. It starts with Binary Coding but
modifies the representation by changing the bit length of the solution
parameters after each iteration. The solution parameters in subsequent
iterations are denoted with respect to the $S_{best}$ from the previous
iteration. The distance from $S_{best}$ is called the Delta Value. This
creates a new hypercube with $S_{best}$ as its origin at each iteration which
is then searched using the crossover operators. In order to maintain diversity
in the population in order to explore the search space, we re-initialize the
population in each iteration. See Algorithm 6 [Mathias and Whitley (1994)].

1 Initialize the population randomly $P(0)$

2 while  _trials $<$ max.Trials $\&$ fitness $>$ threshold_ do

3 while  _Diversity Metric $>$ 1_ do

4 Apply Crossover Operators to get candidates for $P(t+1)$

5 Evaluate fitness function for each individual

6 Insert valid off-springs to $P(t+1)$

7

8 end while

9 Set the best new solution as $S_{best}$

10 Re-initialize Population

11 if  _Delta Value $==0$_ then

12 if  _Parameter Length $>$ Original Parameter Length_ then

13 Encode parameter using $1$ additional bit

14 end if

15

16 else

17 if  _Parameter Length $>$ Lower Bound_ then

18 Encode parameter using $1$ less bit

19 end if

20

21 end if

22

23 end while

Algorithm 6 Delta Coding

####  Fuzzy Coding

In fuzzy coding, for each parameter of the solution we choose one or more from
a set of fuzzy sets (NM, Ms, ZR, PS, PM). Each solution parameter has a
membership function (fuzzy sigmoid, fuzzy gaussian or fuzzy normal). These
membership functions are used to choose values from selected fuzzy sets
corresponding to the solution parameter. The value of the solution parameter
is denoted by the weighted mean of these chosen vales as per the membership
function weights. The membership functions for the fuzzy set NM are given as
follows [Sharma and Irwin (2003)]:

$\displaystyle\text{Fuzzy Sigmoid}:y$
$\displaystyle=-ae^{-(\frac{x-c}{\frac{d}{3}})^{2}}-b-1$
$\displaystyle\text{Fuzzy Gaussian}:y$
$\displaystyle=ae^{-(\frac{x-c}{\frac{d}{3}})^{2}}+b$
$\displaystyle\text{Fuzzy Normal}:y$ $\displaystyle=1-\frac{x-c}{d}$

IN this coding technique, the genetic algorithm need only optimize the fuzzy
sets and the membership function corresponding to each solution parameter. As
the actual value of the parameter is indirectly obtained using the chosen
fuzzy sets and the membership function, the genetic operators (crossover and
mutation) are independent of the actual values. We only need to specify the
type of the parameter and its range beforehand.

###  Selection and Replacement

In every iteration, we need to choose the set of parents from the population,
this can be done by using various policies. Some of the important selection
policies are Roulette Wheel Selection, Ranking Selection, Tournament Selection
and Elitist Selection.

####  Roulette Wheel Selection

In this policy, the selection probability of an individual is proportional to
its fitness value. We can think of this policy as a roulette wheel which is
partitioned as per the fitness of each individual.

$\displaystyle P_{i}$
$\displaystyle=\frac{f(x_{i})}{\sum_{i=1}^{N_{P}}f(x_{i})}$ $\displaystyle i$
$\displaystyle=1,2,...,N_{P}$

####  Ranking Selection

In this policy, the individuals are sorted in ascending order according to
their fitness and the best one is assigned rank 1. This policy applies a
uniform selective pressure on all the individuals unlike the roulette wheel
policy. The probability of the selection of $i^{th}$ individual as a parent is
given by:

$\displaystyle P_{i}$
$\displaystyle=\frac{1}{N_{P}}(\beta-2(\beta-1)\frac{i-1}{N_{P}-1})$
$\displaystyle\beta$ $\displaystyle\in[0,2]$

####  Tournament Selection

In this policy two individuals $i$ and $j$ are picked up randomly with
replacement. The probability that the individual $i$ is selected is given by:

$\displaystyle P_{i}=\frac{1}{1+exp(\frac{f_{j}-f_{i}}{T})}$

We stop after the required number of parents have been selected.

####  Replacement

After choosing the parents and obtaining the offspring, we replace some
individuals from the population. There are several policies to accomplish this
too - complete generational replacement, random replacement, replace worst
policy, replace oldest policy and replce the similar parent policy.

###  Reproduction

Reproduction is simulated by Crossover and Mutation operators. Crossover
operator takes two individuals as parents and gives two new individuals by
swapping some parts of the chromosomes of the parents. Mutation is a unary
operator,it needs only one individual as a parent to create an offspring.

There are different types of crossover operators - one point crossover, two
point crossover, multi-point crossover and uniform crossover. [See Figure 4]
Uniform Crossover swaps are more disruptive and exploratory hence, suitable
for small populations. Two point crossover is generally good for large
populations, It has been demonstrated in many experiments that Two Point
crossover works better than One Point crossover in most of the cases.

Mutation is another genetic operator it randomly chooses one or more bits from
the chromosome and flips them to get a new individual. Mutations introduce the
necessary noise to perform hill climbing in genetic algorithms.

Figure 4: Crossover Operators - a) One Point b) Two Point c) Multi Point d)
Uniform Crossover. Image taken from [Du and Swamy (2016)]

###  Applications of Evolutionary Algorithms

Evolutionary Algorithms have many applications in search and optimization.
Some variants of EAs are clubbed with machine learning algorithms to find
their optimal tuning parameters so as to find the best model. EAs can perform
symbolic regression which are often very difficult for traditional regression
algorithms.

####  Symbolic Regression with Genetic Algorithms

Genetic Programming is a variant of Evolutionary Algorithm which can be used
for symbolic regression. The encoding used in genetic programming is in form
of hierarchical trees. We can represent mathematical function with the help of
these trees. See Figure 5. Now, all weed to perform symbolic regression is to
define the crossover and mutation operators. See Figure 6.

Figure 5: Tree Encoding for $f(x)=x^{2}+sin(x/3)$. Image taken from [Du and
Swamy (2016)] Figure 6: Crossover and Mutation in Symbolic Regressio. Image
taken from [Du and Swamy (2016)]

####  Support Vector Regression with Genetic Algorithms

Chen and Wang (2007) have used Genetic Algorithms for tuning the parameters of
the traditional Support Vector Regression algorithm to forecast tourism data.
Neural Networks have been shown to work very well on this problem. However,
they have some drawbacks such as large number of parameters to train, risk of
overfitting and difficulty in obtaining a stable solution. Support Vector
Regression is also known to work excellently for non-linear regression
problems but it needs careful setting of its parameters to avoid bad
performance and overfitting. But we don’t have any general rule to tune these
parameters, it comes down to knowledge of the experimenter and trial and error
methods for this.

This is the part where genetic algorithms come into picture, Chen and Wang
(2007) applied GA to find the optimal parameters for SVR. Their results are at
par with the Neural Network framework used for the tourism forecasting.

##  Conclusion

We have discussed some very beautiful algorithms inspired by nature which can
be used to solve search and optimization problems. Also, the problems we
discussed are of great practical importance, take Aircraft Conflict Resolution
as an example and Job Shop Scheduling. Solving these problems can save lives
on one hand and time and money on the other. Yet, this is not all, there is a
whole world of such algorithms, algorithms which were designed by observing
and modelling the miracles that nature does. Algorithms like Bee Hive
Optimization, Cuckoo Search, Bacterial Foraging, Glow-worm Swarm Optimization,
Biomolecular Computing and so on, the list is almost endless. All of them
helping us to solve complex problems, mostly problems which are NP Hard.

## References

  * Bianchi et al. (2009) Leonora Bianchi, Marco Dorigo, Luca Maria Gambardella, and Walter J Gutjahr.  A survey on metaheuristics for stochastic combinatorial optimization.  _Natural Computing_ , 8(2):239–287, 2009.
  * Blum and Sampels (2002) Christian Blum and Michael Sampels.  Ant colony optimization for fop shop scheduling: a case study on different pheromone representations.  In _Evolutionary Computation, 2002. CEC’02. Proceedings of the 2002 Congress on_ , volume 2, pages 1558–1563. IEEE, 2002.
  * Chen and Wang (2007) Kuan-Yu Chen and Cheng-Hua Wang.  Support vector regression with genetic algorithms in forecasting tourism demand.  _Tourism Management_ , 28(1):215–226, 2007.
  * Clerc (2010) Maurice Clerc.  _Particle swarm optimization_ , volume 93.  John Wiley & Sons, 2010.
  * Dorigo and Blum (2005) Marco Dorigo and Christian Blum.  Ant colony optimization theory: A survey.  _Theoretical computer science_ , 344(2-3):243–278, 2005.
  * Dorigo et al. (1999) Marco Dorigo, Gianni Di Caro, and Luca M Gambardella.  Ant algorithms for discrete optimization.  _Artificial life_ , 5(2):137–172, 1999.
  * Dorigo et al. (2006) Marco Dorigo, Marco Birattari, and Thomas Stutzle.  Artificial ants as a computational intelligence technique.  _IEEE computational intelligence magazine_ , 1:28–39, 2006\.
  * Du and Swamy (2016) Ke-Lin Du and MNS Swamy.  _Search and Optimization by Metaheuristics_.  Springer, 2016.
  * Durand and Alliot (2009) Nicolas Durand and Jean-Marc Alliot.  Ant colony optimization for air traffic conflict resolution.  In _ATM Seminar 2009, 8th USA/Europe Air Traffic Management Research and Developpment Seminar_ , 2009.
  * Handl et al. (2006) Julia Handl, Joshua Knowles, and Marco Dorigo.  Ant-based clustering and topographic mapping.  _Artificial life_ , 12(1):35–62, 2006.
  * Kirkpatrick et al. (1983) Scott Kirkpatrick, C Daniel Gelatt, Mario P Vecchi, et al.  Optimization by simulated annealing.  _science_ , 220(4598):671–680, 1983.
  * Mathias and Whitley (1994) Keith E. Mathias and L. Darrell Whitley.  Changing representations during search: A comparative study of delta coding.  _Evol. Comput._ , 2(3):249–278, September 1994\.  ISSN 1063-6560.  doi: 10.1162/evco.1994.2.3.249.  URL http://dx.doi.org/10.1162/evco.1994.2.3.249.
  * Pezzella et al. (2008) F Pezzella, G Morganti, and G Ciaschetti.  A genetic algorithm for the flexible job-shop scheduling problem.  _Computers & Operations Research_, 35(10):3202–3212, 2008.
  * Poli et al. (2007) Riccardo Poli, James Kennedy, and Tim Blackwell.  Particle swarm optimization.  _Swarm intelligence_ , 1(1):33–57, 2007.
  * Shao et al. (2013) Xinyu Shao, Weiqi Liu, Qiong Liu, and Chaoyong Zhang.  Hybrid discrete particle swarm optimization for multi-objective flexible job-shop scheduling problem.  _International Journal of Advanced Manufacturing Technology_ , 67, 2013.
  * Sharma and Irwin (2003) Sanjay K. Sharma and George W Irwin.  Fuzzy coding of genetic algorithms.  _IEEE Transactions on Evolutionary Computation_ , 7(4):344–355, 2003.
  * Van Laarhoven et al. (1992) Peter JM Van Laarhoven, Emile HL Aarts, and Jan Karel Lenstra.  Job shop scheduling by simulated annealing.  _Operations research_ , 40(1):113–125, 1992.
  * Yamada and Nakano (1997) Takeshi Yamada and Ryohei Nakano.  Job shop scheduling.  _IEE control Engineering series_ , pages 134–134, 1997.