Phase Transition in Ant Colony Optimization

1. Materials and Methods

There are N two-choice quizzes, each of which is answered by a large number of ants sequentially[21]. These quizzes are labeled by $n=1,\cdots ,N$. The answer provided by the t-th ant is denoted as $X(n,t)\in \{0,1\}$, where $X(n,t)=1$ indicates a correct answer, and $X(n,t)=0$ indicates an incorrect answer. Each ant receives 1 point for a correct answer. After ant t has answered all N questions, the total points (TP) earned by the ant can be calculated using the following equation,

$$\mathrm{TP}\left(t\right)=\sum _{n=1}^{N}X(n,t).$$

Ant t deposits pheromones on his answer $X(n,t)\in \{0,1\}$. The amount of the pheromones is $\mathrm{TP}\left(t\right)$. We assume that the pheromones evaporate and decrease by $e^{-1/\tau }$ per unit time. Here $\tau$ represents the time scale of the pheromone evaporation.

Ant $t+1$ observes the total values of the pheromones associated with question m for each choice $X(m,t+1)=\{0,1\}$. We denote the total value of pheromones that remains on all the questions after ant t has answered,

$$S\left(t\right)\equiv \sum _{s=1}^t\mathrm{TP}\left(s\right)e^{t-s}=\sum _{s=1}^t\sum _{n=1}^{N}X(n,s)e^{t-s}.$$

(1)

Then, the remaining pheromone on $X(m,s)=1,1\le s\le t$ is,

$$S(m,t)\equiv \sum _{s=1}^t\mathrm{TP}\left(s\right)X(m,s)e^{t-s}=\sum _{s=1}^t\sum _{n=1}^{N}X(n,s)X(m,s)e^{t-s}.$$

(2)

The remaining pheromone on $X(m,s)=0,1\le s\le t$ is given by $S\left(t\right)-S(m,t)$.

Ants are not highly intelligent, and the probability of them making the correct choice in the two-choice quizzes by themselves is 1/2. The information provided by $TP\left(s\right)$ gives them an indirect clue about the correct choice. If TP$\left(s\right)>N/2$, the posterior probability for $X(m,s)=1$ is larger than $1/2$. Similarly, if $S(m,t)$ is greater than $S\left(t\right)/2$, the posterior probability for $X(m,t+1)=1$ is greater than $1/2$. In ACO, a decision function is introduced that uses the values of the pheromones as follows,

$$P(X(m,t+1)=1)=(1-ϵ)\frac{S{(m,t)}^{\alpha }}{S{(m,t)}^{\alpha }+{(S\left(t\right)-S(m,t))}^{\alpha }}+\frac{1}{2}ϵ.$$

(3)

Here, the exponent $\alpha$ determines the response of the choice to the values of the pheromones. $ϵ>0$ is a small positive constant to avoid the absorbing states $S(m,t)=0$ and $S(m,t)=S\left(t\right)$ of the stochastic process.

We denote the ratio of the remaining pheromones on the correct choices as $Z(m,t)$,

$$Z(m,t)\equiv \frac{S(m,t)}{S\left(t\right)}.$$

We divide both the denominator and numerator of eq.(3) by $S{\left(t\right)}^{\alpha }$, and the probability of the correct choice $X(m,t+1)=1$ is expressed as,

$$P(X(m,t+1)=1)=(1-ϵ)\left(\frac{Z{(m,t)}^{\alpha }}{Z{(m,t)}^{\alpha }+{(1-Z(m,t))}^{\alpha }}\right)+\frac{1}{2}ϵ=f\left(Z(m,t)\right).$$

Here, $f\left(z\right)$ is defined as

$$f\left(z\right)\equiv (1-ϵ)\left(\frac{z^{\alpha }}{z^{\alpha }+{(1-z)}^{\alpha }}\right)+\frac{1}{2}ϵ.$$

(4)

We note that $f(1/2)=1/2$ and $f(1-x)=1-f\left(x\right)$. The slope of $f\left(x\right)$ at $x=1/2$ is $(1-ϵ)\alpha$. We also introduce the discount factor $D\left(t\right)$ and the ratio of correct answers $Z\left(t\right)$ as,

$$\begin{array}{ccc}\hfill D\left(t\right)& =& \sum _{n=1}^N\sum _{s=1}^t{e}^{-(t-s)/\tau }=N\left(\frac{1-e^{-t/\tau }}{1-e^{-1/\tau }}\right)\hfill \\ \hfill Z\left(t\right)& =& \frac{S\left(t\right)}{D\left(t\right)}.\hfill \end{array}$$

1.1. Stochastic Differential equation

First, we derive the recursive relation for $S\left(t\right)$ and $D\left(t\right)$. According to the definition, $D(t+1)$ and $S(t+1)$ obey the next recursive relations,

$$\begin{array}{ccc}\hfill D(t+1)& =& N+D\left(t\right)e^{-1/\tau }\hfill \\ \hfill S(t+1)& =& \sum _{i=1}^{N}X(i,t+1)+S\left(t\right)e^{-1/\tau }.\hfill \end{array}$$

If $\tau$ is finite, for $t>>\tau >>1$, we have $D\left(t\right)\simeq N\tau$. In the limit $\tau \to \infty$, the pheromone does not evaporate, and we have $D\left(t\right)=Nt$.

For $N>>1$, we can replace the sum of $X(i,s),i=1,\cdots ,N$ with the sum of the expected values using the law of large numbers. We obtain,

$$S(t+1)\simeq \sum _{i=1}^{N}f\left(Z(i,t)\right)+S\left(t\right)e^{-1/\tau }.$$

We denote the average of $f\left(Z\right(i,t\left)\right)$ as $\overline{f\left(Z\right(i,t\left)\right)}$.

$$\overline{f\left(Z\right(i,t\left)\right)}\equiv \frac{1}{N}\sum _{i=1}^{N}f\left(Z(i,t)\right).$$

We have

$$S(t+1)\simeq S\left(t\right)e^{-1/\tau }+N\overline{f\left(Z\right(i,t\left)\right)}.$$

The recursive relation for $Z\left(t\right)$ is

$$Z(t+1)\simeq \frac{S(t+1)}{D(t+1)}=\frac{D(t+1)-N}{D(t+1)}·Z\left(t\right)+\frac{N}{D(t+1)}·\overline{f\left(Z\right(i,t)}$$

$\Delta Z\left(t\right)\equiv Z(t+1)-Z\left(t\right)$ is given as,

$$\Delta Z\left(t\right)\simeq \frac{N}{D(t+1)}(\overline{f\left(Z\right(i,t\left)\right)}-Z\left(t\right)).$$

In the continuous time limit $dt=1\to 0$, we obtain

$$\frac{d}{dt}Z\left(t\right)=\frac{N}{D(t+1)}(\overline{f\left(Z\right(i,t\left)\right)}-Z\left(t\right)).$$

(5)

One see that $Z\left(t\right)$ converges to $\overline{f\left(Z\right(i,t\left)\right)}$. However, when the pheromones evaporate and $\tau <\infty$, $D\left(t\right)\simeq N\tau$ and the prefactor of the differential equation is $1/\tau$. If one assume that the dynamics of $Z(i,t)$ is faster than that of $Z\left(t\right)$(adiabatic approximation), the time scale of the convergence is given by $\tau$ as $\overline{f\left(Z\right(i,t\left)\right)}-Z\left(t\right)\propto e^{-t/\tau }$. When the pheromone does not evaporate and $\tau \to \infty$, $D\left(t\right)\simeq Nt$. The prefactor of the differential equation is $1/t$, and $\overline{f\left(Z\right(i,t\left)\right)}-Z\left(t\right)\propto t^{-1}$. The convergence becomes extremely slow.

Next, we study the dynamics of $Z(m,t)$. The recursive relation for $S(m,t)$ is,

$$S(m,t+1)=S(m,t)e^{-1/\tau }+X(m,t+1)(\sum _{i=1\ne m}^{N}X(i,t+1)+1).$$

$Z(m,t+1)$ is then estimated as,

$$\begin{array}{ccc}\hfill Z(m,t+1)& =& \frac{S(m,t+1)}{S(t+1)}=\frac{S\left(t\right)e^{-1/\tau }}{S(t+1)}·Z(m,t)+\frac{{\sum }_{i\ne m}X(i,t+1)+1}{S(t+1)}·X(m,t+1)\hfill \\ & \simeq & \frac{S(t+1)-N\overline{f\left(Z\right(i,t\left)\right)}}{S(t+1)}·Z(m,t)+\frac{N\overline{f\left(Z\right(i,t\left)\right)}+(1-f\left(Z(m,t)\right)}{S(t+1)}·X(m,t+1).\hfill \end{array}$$

Using $S(t+1)=D(t+1)Z(t+1)$, we obtain

$$\Delta Z(m,t)=\frac{N\overline{f\left(Z\right(i,t\left)\right)}}{D(t+1)Z(t+1)}\left(\left(1+\frac{(1-f(Z(m,t)\left)\right)}{N\overline{f\left(Z\right(i,t\left)\right)}}\right)X(m,t+1)-Z(m,t)\right).$$

We denote the history of $Z\left(s\right),{\left\{Z(i,s)\right\}}_{i=1,\cdots ,N},s=1,\cdots ,t$ as $H_t$, and the conditional expected value of $\Delta Z(m,t)$ is estimated as

$$E(\Delta Z(m,t)|H_t)=\frac{N\overline{f\left(Z\right(i,t\left)\right)}}{D(t+1)Z(t+1)}\left(\left(1+\frac{(1-f(Z(m,t)\left)\right)}{N\overline{f\left(Z\right(i,t\left)\right)}}\right)f\left(Z(m,t)\right)-Z(m,t)\right).$$

(6)

Likewise, the conditional variance of $\Delta Z(m,t)$ can be approximated as,

$$V(\Delta Z(m,t)|H_t)\simeq {\left(\frac{N\overline{f\left(Z\right(i,t)}}{D(t+1)Z(t+1)}\right)}^{2}f\left(Z(m,t)\right)(1-f\left(Z(m,t)\right)).$$

Here, we neglect the subleading terms in ${\left(1+\frac{(1-f(Z(m,t)\left)\right)}{N\overline{f\left(Z\right(i,t\left)\right)}}\right)}^2$. We read the drift and diffusion term from the results and the SDEs are,

$$dZ(m,t)=E(\Delta Z(m,t)|H_t)dt+\sqrt{V(\Delta Z(m,t)|H_t)}dW\left(t\right),m=1,\cdots ,N.$$

(7)

Here, $W\left(t\right)$ is the Wiener process. Eq.(5) and eq.(7) describe the dynamics of the system. The system can be described as a multi-variate Pólya urn process.

We note that $\frac{(1-f(Z(m,t)\left)\right)}{N\overline{f\left(Z\right(i,t\left)\right)}}$ in eq.(6) breaks the $Z_2$ symmetry of the system. If one neglects the term, $E(\Delta Z(m,t)|H_t)$ is proportional to $f\left(Z\right(m,t\left)\right)-Z(m,t)$. As $f(1-x)=1-f\left(x\right)$, $f(1-x)-(1-x)=-\left(f\right(x)-x)$ holds. $Z(m,t)$ and $1-Z(m,t)$ obeys the same dynamics and we call the symmetry $Z_2$ symmetry. The term is always positive and drives $Z(m,t)$ in the positive direction. As the term is proportional to $1/N$, the strength of the $Z_2$-symmetry breaking field becomes smaller as N becomes larger.

2. Results

We analyze the SDEs given in eq.(7) and investigate the convergence properties of $Z(m,t)$. As the convergence behavior relies on the initial value of $Z(m,t=t_0)$ in the context of the non-linear Pólya urn model, we commence by examining the distribution of $Z(m,t)$.

2.1. Initial distribution of $\left\{Z\right(m,t\left)\right\}$

We assume that ants adopt $\alpha =0$ and do not respond to the values of the pheromones for $t=1,\cdots ,t_0$. The ants answer the questions independently and $P\left(X\right(i,t)=1)=1/2$. We estimate $Z\left(t\right)$ and $Z(m,t)$ for $\tau <t\le t_0$ as follows.

$$\begin{array}{ccc}\hfill Z\left(t\right)& \sim & \mathrm{N}\left(\frac{1}{2},\frac{D_h\left(t\right)}{4D{\left(t\right)}^2}\right),\hfill \\ \hfill Z(m,t)& \sim & \mathrm{N}\left(\frac{1}{2}+\frac{1}{2N},\frac{N{D}_h\left(t\right)}{4D{\left(t\right)}^2}\right),\hfill \\ \hfill D_h\left(t\right)& \equiv & N\sum _{s=1}^t{e}^{-2(t-s)/\tau }=N·\frac{1-e^{-2t/\tau }}{1-e^{-2/\tau }}.\hfill \end{array}$$

(8)

The details of the calculations are given in Appendix A. If $\tau$ is finite, we have $D_h\left(t\right)\simeq N\tau /2$ for $t>>\tau >>1$. In the limit $\tau \to \infty$, the pheromone does not evaporate and we have $D_h\left(t\right)=Nt$.

The essential differences between $Z\left(t\right)$ and $Z(m,t)$ include a shift of the expected value by $1/2N$ and the presence of a factor of N in the numerator of the variance of $Z(m,t)$. The shift of $1/2N$ arises from the fact that $X(m,s)$ in $S(m,t)$ is 1, which is larger than $E\left[X\right(i,s\left)\right]=1/2$ for $i\ne m$. The value of the pheromone contains information about the correct choice, leading to $E\left[S(m,t)\right]>\frac{1}{2}E\left[S\left(t\right)\right]$. However, in the "cheating" process, the variables $X(i,s),i=1,\cdots ,N$ are combined by $X(m,s)$ as in eq.(2), resulting in a larger variance for $S(m,t)$. The factor of N in the numerator of the variance of $Z(m,t)$ is a consequence of this combination process.

From the distribution of $Z(m,t_0)$, one can determine the values of $\tau$ or $t_0$ that guarantee that $Z(m,t_0)$ is greater than 0.5 for the limits $t>>\tau >>1$ and $\tau \to \infty$, respectively. With a confidence level of 1%, $\tau$ and $t_0$ should satisfy the following conditions:

$$\begin{array}{ccc}& & t>>\tau >>1:\frac{1}{2N}\ge \frac{2.58}{2\sqrt{2}\tau }\to \tau \ge 3.33N^2\hfill \\ & & \tau \to \infty :\frac{1}{2N}\ge \frac{2.58}{2\sqrt{t_0}}\to t_0\ge 6.66N^2.\hfill \end{array}$$

2.2. $\tau \to \infty$ case

We take the limit $\tau \to \infty$ in eq.(5) and eq.(7). We replace $D(t+1)$ and $D_h(t+1)$ with $N(t+1)$ and $N(t+1)$, respectively. This results in the following equations:

$$\begin{array}{ccc}\hfill dZ(m,t)& =& \frac{\overline{f\left(Z\right(i,t\left)\right)}}{(t+1)Z(t+1)}\left(\left(1+\frac{(1-f(Z(m,t)\left)\right)}{N\overline{f\left(Z\right(i,t\left)\right)}}\right)f\left(Z(m,t)\right)-Z(m,t)\right)dt\hfill \\ & +& \left(\frac{\overline{f\left(Z\right(i,t\left)\right)}}{(t+1)Z(t+1)}\right)\sqrt{f\left(Z\right(m,t\left)\right(1-f\left(Z\right(m,t\left)\right))}dW\left(t\right),m=1,\cdots ,N\hfill \\ \hfill \frac{d}{dt}Z\left(t\right)& =& \frac{1}{(t+1)}(\overline{f\left(Z\right(i,t\left)\right)}-Z\left(t\right)).\hfill \end{array}$$

The initial conditions of $Z\left(t\right)$ and $Z(m,t)$ at $t=t_0$ are as follows:

$$\begin{array}{ccc}\hfill Z\left(t_0\right)& \sim & \mathrm{N}\left(\frac{1}{2},\frac{1}{4N{t}_0}\right),\hfill \\ \hfill Z(m,t_0)& \sim & \mathrm{N}\left(\frac{1}{2}+\frac{1}{2N},\frac{1}{4t_0}\right).\hfill \end{array}$$

In the adiabatic approximation, where the time development of $Z(m,t)$ is much faster than that of $Z\left(t\right)$, the structure of the SDE for each $Z(m,t)$, where $m=1,\cdots ,N$, is the same as that of a non-linear P’olya urn [22,23]. The probability for $Z(m,t)$ to converge to a stable solution of the following equation is positive[24].

$$\left(1+\frac{(1-f(z\left)\right)}{N\overline{f\left(Z\right(i,t\left)\right)}}\right)f\left(z\right)=z.$$

(9)

Here, a stable (unstable) solution of eq.(9) means that the curve of the left-hand-side of the equation crosses the diagonal curve $y=z$ in in the downward (upward) direction[24]. We denote the stable and unstable solutions as $z_s$ and $z_u$, respectively.

In the limit $N\to \infty$, the positive driving force $\frac{(1-f(z\left)\right)}{N\overline{f\left(Z\right(i,t\left)\right)}}$ in eq.(9) disappears. The system becomes $Z_2$ symmetric and $z=1/2$ is a the solution as $f(1/2)=1/2$. If N is finite, the positive driving force breaks the $Z_2$ symmetry. We plot $(1+b(1-f\left(z\right)\left)f\right(z)$ in Figure 1. $b>0$ corresponds with $1/N\overline{f\left(Z\right(i,t\left)\right)}$ and $b\simeq 0.01\sim 0.02$ for $N=100$ and $\overline{f\left(Z\right(i,t)}=1/2\sim 1$. In the left and the right figure, we adopt $b=0,ϵ=0.1$ and $b=0.2,ϵ=0.1$, respectively.

In the $Z_2$ symmetric case ($b=0$), the stability of the solution $z=1/2$ depends on the slope of $f\left(z\right)$ at $z=1/2$. As $f^{\prime }(1/2)=(1-ϵ)\alpha$, the critical value of $\alpha$ is ${\alpha }_c=1/(1-ϵ)$. If $\alpha <(>){\alpha }_c$, $z=1/2$ is (un)stable. If $\alpha ={\alpha }_c$, the curve of $f\left(z\right)$ is tangential to the diagonal. $Z(m,t)$ converges to $1/2$ for $\alpha \le {\alpha }_c$, as $z=1/2$ is the unique stable solution $z_s$. For $\alpha >{\alpha }_c$, there appears two stable solutions $z_s,z_s^{\prime }$, one($z_s^{\prime }$) is in $(0,0.5)$ and the other($z_s$) is in $(0.5,1.0)$. $z=0.5$ becomes the unstale solution $z_u$.

The dotted lines Figure 2 shows the solutions vs. $\alpha$ for the $Z_2$ symmetric case. For $\alpha <{\alpha }_c$, $z=1/2$(black dotted line) is the stable solution. For $\alpha >\alpha >c$, $z=1/2$ (gray dotted line) becomes unstable ($z_u=1/2$) and two stable solution $z_s,z_s^{\prime }$ departs from $z=1/2$ continuously with $\alpha >{\alpha }_c$. Which stable solution does $Z(m,t)$ converge depends on the initial value of $Z(m,t_0)$. In general, if $Z(m,t_0)$ is greater (smaller) than $1/2$, the probability of the convergence to 1 is greater (smaller) than $1/2$. $z_u$ determines the "attractive domains" for the stable solutions $z_s,z_s^{\prime }$. The susceptibility of the expected value of $Z(m,t)$ to the initial value $Z(m,t_0)$ is the order parameter of the non-linear Pólya urn[22,23]. As the order parameter is proportional to the difference of the two stable states, the order parameter is a continuous function of $\alpha$ and the phase transition is continuous.

In the $Z_2$ asymmetric case ($b\ne 0$), for small values of $\alpha$, there is a stable solution $z_s$ in the range $(0.5,1)$. $z_s$ increase with $\alpha$ and at some critical value ${\alpha }_c$ of $\alpha$, the curve $f\left(z\right)(1+b(1-f\left(z\right))$ becomes tangential to the diagonal at $z=z_t$. $z_t$ is known as touchpoint and to be stable[25]. As $\alpha$ continues to increase beyond ${\alpha }_c$, two significant changes occur: a new stable solution $z_s^{\prime }$ emerges from the touchpoint $z_t$, while an unstable solution $z_u$ also becomes apparent.

When $\alpha <{\alpha }_c$, only one stable solution, $z_s$, exists within the range $1/2<z_s<1$. Conversely, for $\alpha >{\alpha }_c$, the specific stable solution to which $Z(m,t)$ converges depends on the initial values of $Z(m,t_0)$. At $\alpha ={\alpha }_c$, both the stable fixed point $z_s$ and the touchpoint $z_t$ remain stable. Which solution $Z(m,t)$ converges to is determined by the initial values of $Z(m,t_0)$ in this case as well. For $\alpha >{\alpha }_c$, the situation mirrors that of $\alpha ={\alpha }_c$. Once $\alpha \ge {\alpha }_c$, the order parameter turns positive, and the phase transition becomes discontinuous.

Figure 3 shows the results of the numerical studies in the limit $\tau \to \infty$. We sampled a trajectory of $Z(m,t)$ and $Z\left(t\right)$ for $1\le t\le {10}^9$ with $N={10}^2$ and $ϵ=0.01$. In the left figure, we present the distribution of $Z(m,t_0)$ for two different values of $t_0$, namely, $t_0\in {10}^3,{10}^6$. The mean value of $Z(m,t_0)$ is approximately $1/2+1/2N$, which aligns with the theoretical predictions. The variance of $Z(m,t_0)$ is given by $1/4t_0$, so the variance for $t_0={10}^3$ is about ${10}^3$ times larger than that for $t_0={10}^6$. Consequently, if we choose $t_0={10}^3$, a significant proportion of $m\in 1,\cdots ,N$ will have $Z(m,t_0)<z_u\simeq 0.5$, resulting in a high probability that $Z(m,t)$ converges to $z_s^{\prime }<1/2$ for $\alpha =2.0$. In such cases, $Z\left(t\right)$ cannot reach 1 due to the convergence of $Z(m,t)$ to $z_s^{\prime }$. On the other hand, if we set $t_0={10}^5$, the ratio of $m\in 1,\cdots ,N$ with $Z(m,t_0)<z_u\simeq 0.5$ is zero, ensuring that $Z(m,t)$ always converges to $z_s>1/2$. As a result, $Z\left(t\right)$ monotonically increases towards 1 for $\alpha =2.0$. For $\alpha =1.0<{\alpha }_c$, where only one stable state $z_s\simeq 1$ exists, $Z(m,t)$ consistently converges to $z_s$. It’s evident that $Z\left(t\right)$ monotonically approaches $z_s$ with time for both $t_0={10}^3$ and $t_0={10}^5$ cases within the range $t\le {10}^9$. In the case of $\alpha =0.5$, where $z_s\simeq 0.5$, $Z\left(t\right)$ experiences relatively little change

2.3. $\tau <\infty$ case

In the case where $\tau$ is finite, we make the assumption that $t\gg \tau \gg 1$ and replace $D(t+1)$ with $N\tau$ in eq.(5) and eq.(7). This leads to the following equations:

$$\begin{array}{ccc}\hfill dZ(m,t)& =& \frac{\overline{f\left(Z\right(i,t\left)\right)}}{\tau Z(t+1)}\left(\left(1+\frac{(1-f(Z(m,t)\left)\right)}{N\overline{f\left(Z\right(i,t\left)\right)}}\right)f\left(Z(m,t)\right)-Z(m,t)\right)dt\hfill \\ & +& \left(\frac{\overline{f\left(Z\right(i,t\left)\right)}}{\tau Z(t+1)}\right)\sqrt{f\left(Z\right(m,t\left)\right)(1-f(Z(m,t))}dW\left(t\right),m=1,\cdots ,N\hfill \\ \hfill \frac{d}{dt}Z\left(t\right)& =& \frac{1}{\tau }(\overline{f\left(Z\right(i,t\left)\right)}-Z\left(t\right)).\hfill \end{array}$$

(10)

The initial conditions for $Z\left(t\right)$ and $Z(m,t)$ at $t=t_0\gg \tau$ are given as follows:

$$\begin{array}{ccc}\hfill Z\left(t_0\right)& \sim & \mathrm{N}\left(\frac{1}{2},\frac{1}{8N\tau }\right),\hfill \\ \hfill Z(m,t_0)& \sim & \mathrm{N}\left(\frac{1}{2}+\frac{1}{2N},\frac{1}{8\tau }\right).\hfill \end{array}$$

(11)

The dynamics of ${Z(m,t)}_{m=1,\cdots ,N}$ are coupled through $Z\left(t\right)$ and $\overline{f}\left(Z(i,t)\right)$. To simplify and analyze this coupled system, we focus on the stationary state of $Z(m,t)$ in the limit $t\to \infty$.

We anticipate that $Z(m,t)m=1,\cdots ,N$ will fluctuate around the stable fixed points of $f\left(z\right)$ in the stationary state. As we observed earlier in the case of $\tau \to \infty$, for $\alpha <{\alpha }_c$, there is only one stable fixed point, and for $\alpha \ge {\alpha }_c$, two stable fixed points exist, one of which is near 1. The stationary distribution is unimodal for $\alpha <{\alpha }_c$ and bimodal for $\alpha \ge {\alpha }_c$. We denote the stationary distribution and the mean value of $Z(m,t)$ as $Pst\left(z\right)$ and ${\mu }_{st}$, respectively. As $f\left(Z\right(m,t\left)\right)$ is the probability for $X(m,t+1)=1$, we can assume $\overline{f\left(Z\right(i,t)}=Z\left(t\right)={\mu }_{st}$ in the stationary state. The SDEs in eq.(10) can be simplified as follows when replacing $\overline{f}\left(Z(i,t)\right)$ and $Z\left(t\right)$ with ${\mu }_{st}$:

$$\begin{array}{ccc}\hfill dZ(m,t)& =& \frac{1}{\tau }\left(\left(1+\frac{(1-f(Z(m,t)\left)\right)}{N{\mu }_{st}}\right)f\left(Z(m,1)\right)-Z(m,t)\right)dt\hfill \\ & +& \frac{1}{\tau }\sqrt{f\left(Z\right(m,t\left)\right)(1-f(Z(m,t))}dW\left(t\right)\hfill \\ & =& A\left(Z\right(m,t\left)\right)dt+B\left(Z\right(m,t\left)\right)dW\left(t\right)\hfill \\ \hfill A\left(z\right)& =& \frac{1}{\tau }\left(\left(1+\frac{(1-f(z\left)\right)}{N{\mu }_{st}}\right)f\left(z\right)-z\right)\hfill \\ \hfill B\left(z\right)& =& \frac{1}{\tau }\sqrt{f\left(z\right)(1-f(z)}\hfill \end{array}$$

(12)

The stationary state with reflecting boundary conditions is determined by a potential solution[26], which can be expressed as:

$$\begin{array}{ccc}\hfill P_{st}\left(z\right)& \propto & \frac{1}{B{\left(z\right)}^2}exp\left({\int }_{1/2}^z\frac{2A\left(y\right)}{B{\left(y\right)}^2}dy\right).\hfill \\ \hfill \frac{2A\left(y\right)}{B{\left(y\right)}^2}& =& 2\tau \frac{f\left(y\right)-y}{f\left(y\right)(1-f(y\left)\right)}+\frac{2\tau }{N{\mu }_{st}}.\hfill \end{array}$$

(13)

The second term, $\frac{2\tau }{N{\mu }_{st}}$, arises from the $Z_2$ symmetry-breaking field and causes a shift in the stationary distribution in the positive direction.

Figure 4 shows $P_{st}\left(z\right)$ in eq.(13) for $ϵ=0.01,N={10}^2$. ${\mu }_{st}$ is chosen so that the mean value of $P_{st}\left(z\right)$ coincides with ${\mu }_{st}$.

$${\mu }_{st}={\int }_0^1{P}_{st}\left(z\right)zdz.$$

The parameters $(\alpha ,{\mu }_{st})$ are $(0.0,0.50)$,$(0.5,0.51)$,$(0.9,0.54)$,$(0.99,0.67)$, $(1/0.99,0.74)$ and $(2.0,0.54)$. As $\alpha$ increases, the peak position shifts in the positive direction, which can be expected by the dependence of the stable solution $z_s$ on $\alpha$ in Figure 2. If $\alpha =1/(1-ϵ)$, the peak appears at $z=1$, since there is only one stable fixed point near 1 in Figure 1. When $\alpha =2$, there are two stable fixed point and the stationary distribution is bimodal.

In order to derive the dependence of ${\mu }_{st}$ and the variance of $P_{st}\left(z\right)$ on $\alpha$, we assume that $Z(m,t)$ fluctuates around ${\mu }_{st}\simeq \frac{1}{2}$ for $\alpha <{\alpha }_c$. We linearrize $f\left(z\right)$ in the vicinity of $z=1/2$ as,

$$f\left(z\right)=\frac{1}{2}+(1-ϵ)\alpha \left(z-\frac{1}{2}\right).$$

We also approximate $B{\left(z\right)}^2$ as ${\mu }_{st}(1-{\mu }_{st})=1/4$, $P_{st}\left(z\right)$ becomes

$$P_{st}\left(z\right)\propto exp\left(-\frac{{\left(z-\left(\frac{1}{2}+\frac{1}{2N(1-(1-ϵ\left)\alpha \right)}\right)\right)}^2}{\left(\frac{2}{8\tau (1-(1-ϵ\left)\alpha \right)}\right)}\right)$$

In the case $\alpha =0$, the ants does not observe the information of the pheromones and decide by themselves. The expected value and the variance are consistent with the results for the initial state in eq.(11). The expected value and the variance increase with $\alpha$ for $0\le \alpha <1/(1-ϵ)$.

The shape of the stationary distribution changes from the monomodal shape for $\alpha <{\alpha }_c$ to the bimodal shape for $\alpha \ge {\alpha }_c$. Figure 5 shows the stationary distribution of $Z(m,t)$ for $\alpha \in \{0.0,0.5,0.9,0.99,1/0.99,2.0\},ϵ=0.01,\tau =100$ and $N={10}^2$. We also plot $P_{st}\left(z\right)$ in eq.(13) with solid line curves. Except for the $\alpha =2.0$ case, the numerical results agree with the theoretical ones. As $\alpha$ increases from 0 to $1/0.99$, the mean value and the variance of $Z(m,t)$ increases. For $\alpha =1/0.99$, the distribution of $Z(m,t)$ has a peak at $z=1$. The distribution becomes bimodal and has two peaks near $z=0$ and $z=1$ for $\alpha =2$. For $\alpha =2.0$, $P_{st}\left(z\right)$ becomes bimodal and the equilibriation time to reach the stationary state becomes extremely long. We think this is the reason for the discrepancy between the numerical and theoretical results.

3. Discussion

We have studied a simple model for ACO and the convergence properties of the solutions. Ants answer many two-choice quizzes in sequence and deposit pheromone as they choose. As the amount of the pheromones is the number of correct answers, the following ants can receive information or hints about the correct choices. We have shown that the model reduces to a multi-variate non-linear Pólya urn process and the pheromones break the $Z_2$ symmetry of the process. By varying the exponent $\alpha$ of the decision function of the ants, there occurs a phase transition about the convergence of the probability of choosing the correct answer for each question in the limit $\tau \to \infty$. For $\tau <\infty$, the change of the stationary distribution between the monomodal and the bimodal shape occurs as we vary $\alpha$.

previous studies have adopted values of $\alpha =1$ or smaller in solving real problems, like TSP. In $\alpha$-annealing, $\alpha$ increases gradually, as shown in previous research[17]. In our study, we have shown that the duration of the period $\alpha =0$ should be long enough to ensure that the initial value of $Z(m,t)$ is in the attractive domain of the good stable state ($Z(m,t)>z_u$) in the case $\tau =\infty$. Subsequently, with $\alpha \ge 1$ in effect, $Z\left(t\right)$ converges to a value close to 1. In the case of $\tau <\infty$, the timescale for pheromone evaporation, represented by $\tau$, should be sufficiently long to maintain the same initial conditions. However, $Z(m,t)$ does not converge to a specific value; instead, it follows a stationary distribution that exhibits both bimodal and monomodal shapes depending on the value of $\alpha$. To achieve a distribution of $Z\left(t\right)$ with a prominent peak near $z=1$, the $\alpha$-annealing process is an effective strategy. Both the stable solution $z_s$ and the distribution of $Z(m,t)$ suggest that, after a lengthy period with $\alpha =0$, it is advantageous to gradually increase $\alpha$ from 1. However, the efficiency of the annealing process depends on the specific problem being addressed. Future study should clarify the efficient $\alpha$-annealing schedule.