# QPSO-CD: Quantum-behaved Particle Swarm Optimization Algorithm with Cauchy
Distribution

Amandeep Singh Bhatia${}^{1,3,^{*}}$, Mandeep Kaur Saggi2, Shenggen Zheng3,
Soumya Ranjan Nayak4
1Chitkara University Institute of Engineering & Technology, Chitkara
University, Punjab, India
2Department of Computer Science, Thapar Institute of Engineering & Technology,
India
3Center for Quantum Computing, Peng Cheng Laboratory, Shenzhen, China
4Amity School Of Engineering & Technology, Amity University Uttar Pradesh,
Noida, India
E-mail: amandeepbhatia.singh@gmail.com${}^{1,^{*}}$

###### Abstract

Motivated by the particle swarm optimization (PSO) and quantum computing
theory, we have presented a quantum variant of PSO (QPSO) mutated with Cauchy
operator and natural selection mechanism (QPSO-CD) from evolutionary
computations. The performance of proposed hybrid quantum-behaved particle
swarm optimization with Cauchy distribution (QPSO-CD) is investigated and
compared with its counterparts based on a set of benchmark problems. Moreover,
QPSO-CD is employed in well-studied constrained engineering problems to
investigate its applicability. Further, the correctness and time complexity of
QPSO-CD are analysed and compared with the classical PSO. It has been proved
that QPSO-CD handles such real-life problems efficiently and can attain
superior solutions in most of the problems. The experimental results shown
that QPSO associated with Cauchy distribution and natural selection strategy
outperforms other variants in context of stability and convergence.

Keywords: Quantum-behaved algorithm, Particle swarm optimization, Engineering
design problems, Cauchy distribution, Quantum computing.

##  I Introduction

In the late 19th century, the theory of classical mechanics experienced
several issues in reporting the physical phenomena of light masses and high
velocity microscopic particles. In 1920s, Bohr’s atomic theory bohr1928quantum
, Heisenberg’s discovery of quantum mechanics robertson1929uncertainty  and
Schrödinger’s wessels1979schrodinger  discovery of wave mechanics influence
the conception of a new field i.e. the quantum mechanics. In 1982, Feynman
feynman1982simulating  stated that quantum mechanical systems can be simulated
by quantum computers in exponential time, i.e. better than with classical
computers. Till then, the concept of quantum computing was thought to be only
a theoretical possibility, but over the last three decades the research has
evolved such as to make quantum computing applications a realistic possibility
wang2012handbook .

In the last two decades, the field of swarm intelligence has got overwhelming
response among research communities. It is inspired by nature and aims to
build decentralized and self-organized systems by collective behavior of
individual agents with each other and with their environment. The research
foundation of swarm intelligence is constructed mostly upon two families of
optimization algorithms i.e. ant colony optimization (Dorigo et at.
dorigo1999ant  and Colorni et al. colorni1992distributed  1992); and particle
swarm optimization (PSO) (Kennedy & Eberhart kennedy1995particle  1995).
Originally, the swarm intelligence is inspired by certain natural behaviors of
flocks of birds and swarms of ants.

In the mid 1990s, particle swarm optimization technique was introduced for
continuous optimization, motivated by flocking of birds. The evolution of PSO
based bio-inspired techniques has been in an expedite development in the last
two decades. It has got attention from different fields such as inventory
planning wang2014modified , power systems alrashidi2009survey , manufacturing
yildiz2009novel , communication networks latiff2007energy , support vector
machines lin2008particle , to estimate binary inspiral signal wang2010particle
, gravitational waves normandin2018particle  and many more. Similar to
evolutionary genetic algorithm, it is inspired by simulation of social
behavior, where each individual is called particle, and group of individuals
is called swarm. In multi-dimensional search space, the position and velocity
of each particle represent a probable solution. Particles fly around in a
search space seeking potential solution. At each iteration, each particle
adjusts its position according to the goal of its own and its neighbors. Each
particle in a neighborhood share the information with others sun2004particle .
Later, each particle keeps the record of best solution experienced so far to
update their positions and adjust their velocities accordingly.

Figure 1: Particles movement in PSO and QPSO algorithm

Since the first PSO algorithm proposed, the several PSO algorithms have been
introduced with plethora of alterations. Recently, the combination of quantum
computing, mathematics and computer science have inspired the creation of
optimization techniques. Initially, Narayanan and Moore narayanan1996quantum
introduced quantum-inspired genetic algorithm (QGA) in 1995. Later, Sun et al.
sun2004particle  applied the quantum laws of mechanics to PSO and proposed
quantum-inspired particle swarm optimization (QPSO). It is the commencement of
quantum-behaved optimization algorithms, which has subsequently made a
significant impact on the academic and research communities alike.

Recently, Yuanyuan and Xiyu yuanyuan2018quantum  proposed a quantum
evolutionary algorithm to discover communities in complex social networks. Its
applicability is tested on five real social networks and results are compared
with classical algorithms. It has been proved that PSO lacks convergence on
local optima i.e. it is tough for PSO to come out of the local optimum once it
confines into optimal local region. QPSO with mutation operator (QPSO-MO) is
proposed to enhance the diversity to escape from local optimum in search
liu2005quantum . Protopopescu and Barhen protopopescu2002solving  solved set
of global optimization problems efficiently using quantum algorithms. In
future, the proposed algorithm can be integrated with matrix product state
based quantum classifier for supervised learning 40 ; 44 .

In this paper, we have combined QPSO with Cauchy mutation operator to add long
jump ability for global search and natural selection mechanism for elimination
of particles. The results shown that it has great tendency to overcome the
problem of trapping into local search space. Therefore, the proposed hybrid
QPSO strengthened the local and global search ability and outperformed the
other variants of QPSO and PSO due to fast convergence feature.

The illustration of particles movement in PSO and QPSO algorithm is shown in
Fig 1. The big circle at center denotes the particle with the global position
and other circles are particles. The particles located away from global
position are lagged particles. The blue color arrows signify the directions of
other particles and the big red arrows point towards the side in which it goes
with high probability. During iterations, if the lagged particle is unable to
find better position as compared to present global position in PSO, then their
impact is null on the other particles. But, in QPSO, the lagged particles move
with higher probability in the direction of gbest position. Thus, the
contribution of lagged particles is more to the solution in QPSO in comparison
with PSO algorithm.

The organization of rest of this paper is as follows: Sect. 2 is devoted to
prior work. In Sect. 3, the quantum particle swarm optimization is described.
In Sect. 4, the proposed hybrid QPSO algorithm with Cauchy distribution and
natural selection mechanism is presented. The experimental results are plotted
for a set of benchmark problems and compared with several QPSO variants in
Sect. 5. The correctness and time complexity are analyzed in Section 6. QPSO-
CD is applied to three constrained engineering design problems in Sect. 7.
Finally, Sect. 8 is the conclusion.

##  II Prior Work

Since the quantum-behaved particle swarm optimization was proposed, various
revised variants have been emerged. Initially, Sun et al. sun2004particle
applied the concept of quantum computing to PSO and developed a quantum Delta
potential well model for classical PSO sun2004global . It has been shown that
the convergence and performance of QPSO are superior as compared to classical
PSO. The selection and control of parameters can improve its performance,
which is posed as an open problem. Sun et al. sun2007using  tested the
performance of QPSO on constrained and unconstrained problems. It has been
claimed that QPSO is a promising optimization algorithm, which performs better
than classical PSO algorithms. In 2011, Sun et al. sun2011quantum  proposed
QPSO with Gaussian distribution (GAQPSO) with the local attenuator point and
compared its results with several PSO and QPSO counterparts. It has been
proved that GAQPSO is efficient and stable with superior features in quality
and robustness of solutions.

Further, Coelho dos2010gaussian  applied GQPSO to constrained engineering
problems and shown that the simulation results of GQPSO are much closer to the
perfect solution with small standard deviation. Li et al. li2012improved
presented a cooperative QPSO using Monte Carlo method (CQPSO), where particles
cooperate with each other to enhance the performance of original algorithm. It
is implemented on several representative functions and performed better than
the other QPSO algorithms in context of computational cost and quality of
solutions. Peng et al. introduced peng2013quantum  QPSO with Levy probability
distribution and claimed that there are very less chances to be stuck in local
optimum.

Researchers have applied PSO and QPSO to real-life problems and achieved
optimal solutions as compared to existing algorithms. Ali et al. 90  performed
energy-efficient clustering in mobile ad-hoc networks (MANET) with PSO. The
similar approach can be followed to analyse and execute mobility over MANET
with QPSO-CD 88 . Zhisheng zhisheng2010quantum  used QPSO in economic load
dispatch for power system and proved superior to other existing PSO
optimization algorithms. Sun et al. sun2006qpso  applied QPSO for QoS
multicast routing. Firstly, the QoS multicast routing is converted into
constrained integer problems and then effectively solved by QPSO with loop
deletion task. Further, the performance is investigated on random network
topologies. It has been proved that QPSO is more powerful than PSO and genetic
algorithm. Geis and Middendorf proposed PSO with Helix structure for finding
ribonucleic acid (RNA) secondary structures with same structure and low energy
geis2011particle . The QPSO-CD algorithm can be used with two-way quantum
finite automata to model the RNA secondary structures bhatia2018modeling .
Bagheri et al. bagheri2014financial  applied the QPSO for tuning the
parameters of adaptive network-based fuzzy inference system (ANFIS) for
forecasting the financial prices of future market. Davoodi et al.
davoodi2014hybrid  introduced a hybrid improved QPSO with Neldar Mead simplex
method (IQPSO-NM), where NM method is used for tuning purpose of solutions.
Further, the proposed algorithm is applied to solve load flow problems of
power system and acquired the convergence accurately with efficient search
ability. Omkar omkar2009quantum  proposed QPSO for multi-objective design
problems and results are compared with PSO. Recently, Fatemeh et al.
fatemeh2019shuffled  proposed QPSO with shuffled complex evolution (SP-QPSO)
and its performance is demonstrated using five engineering design problems.
Prithi and Sumathi prithi2020ld2fa  integrated the concept of classical PSO
with deterministic finite automata for transmission of data and intrusion
detection. The proposed algorithm QPSO-CD can be used with quantum
computational models for wireless communication 10 ; 20 ; 30 .

##  III Quantum Particle Swarm Optimization

Before we explain our hybrid QPSO-CD algorithm mutated with Cauchy operator
and natural selection method, it is useful to define the notion of quantum
PSO. We assume that the reader is familiar with the concept of classical PSO;
otherwise, reader can refer to particle swarm optimization algorithm
kennedy2010particle ; shi2001particle . The specific principle of quantum PSO
is given as:

In QPSO, the state of a particle can be represented using wave function
$\psi(x,t)$. The probability density function $|\psi(x,t)|^{2}$ is used to
determine the probability of particle occurring in position x at any time t
sun2004particle ; sun2006qpso . The position of particles is updated according
to equations:

$x_{i,j}(t+1)=p_{i,j}(t)\pm\alpha.|mbest_{i,j}(t)-x_{i,j}(t)|.ln(1/u)$ (1)
$p_{i,j}(t)=(\phi.P_{i,j}(t)+(1-\phi).G_{j}(t)),(1\leq i\leq N,1\leq j\leq M)$
(2)

where each particle must converge to its local attractor
$p=(p_{1},p_{2},...,p_{D})$, where D is the dimension, N and M are the number
of particles and iterations respectively, $P_{i,j}$ and $G_{j}$ denote the
previous and optimal position vector of each particle respectively,
$\phi=c_{1}.r_{1}/(c_{1}r_{1}+c_{2}r_{2})$, where $c_{1}$; $c_{2}$ are the
acceleration coefficients, $r_{1}$; $r_{2}$ and u are normally distributed
random numbers in (0, 1), $\alpha$ is contraction-expansion coefficient and
mbest defines the mean of best positions of particles as:

$mbest_{i,j}(t)=\dfrac{1}{N}\sum_{i=1,j=1}^{N,M}P_{i,j}(t)=\bigg{(}\dfrac{1}{N}\sum_{i=1}^{N}P_{i,1}(t),\dfrac{1}{N}\sum_{i=1}^{N}P_{i,2}(t),...,\dfrac{1}{N}\sum_{i=1}^{N}P_{i,D}(t)\bigg{)}$
(3)

In Eq. (1), $\alpha$ denotes contraction-expansion coefficient, which is setup
manually to control the speed of convergence. It can be decreased linearly or
fixed. In PSO, $\alpha<1.782$ to ensure convergence performance of the
particle. In QPSO-CD, the value of $\alpha$ is determined by
$\alpha$=1-(1.0-0.5)k/M, i.e. decreases linearly from 1.0 to 0.5 to attain
good performance, where k is present iteration and M is maximum number of
iterations.

##  IV Hybrid Particle Swarm Optimization

The hybrid quantum-behaved PSO algorithm with Cauchy distribution and natural
selection strategy (QPSO-CD) is described as follows:

The QPSO-CD algorithm begins with the standard QPSO using equations (1), (2)
and (3). The position and velocity of particles cannot be determined exactly
due to varying dynamic behavior. So, it can only be learned with the
probability density function. Each particle can be mutated with Gaussian or
Cauchy distribution. We mutated QPSO with Cauchy operator due to its ability
to make larger perturbation. Therefore, there is a higher probability with
Cauchy as compared to Gaussian distribution to come out of the local optima
region. The QPSO algorithm is mutated with Cauchy distribution to increase its
diversity, where mbest or global best position is mutated with fixed mutation
probability (Pr). The probability density function $(f(x))$ of the standard
Cauchy distribution is given as:

$f(x)=\dfrac{1}{\pi(1+x^{2})}~{}~{}~{}~{}-\infty<x<\infty,$ (4)

It should be noted that mutation operation is executed on each vector by
adding Cauchy distribution random value (D(.)) independently such that

$x^{{}^{\prime}}=x+\phi D(.)$ (5)

where $x^{{}^{\prime}}$ is new location after mutated with random value to x.
At last, the position of particle is selected and the particles of swarm are
sorted on the basis of their fitness values after each iteration. Further,
substitute the group of particles having worst fitness values with the best
ones and optimal solution is determined. The main objective of using natural
mechanism is to refine the capability and accuracy of QPSO algorithm.

Algorithm 1 QPSO-CD algorithm

1:The swarm is initialized with random numbers distributed uniformly: random
$x_{i}$.

2:Do

3:$\alpha$ decreases linearly from 1.0 to 0.5

4:for k=1 to M do do

5: $\alpha=1-(1.0-0.5).k/M$

6: if $Pr<rand(0,1)$ then

7: Calculate the mbest of the swarm using Eq (3)

8: end if

9: for i=1: to N do do

10: if Fitness$(p_{i})<$ Fitness$(x_{i})$ then $p_{i}=x_{i}$;

11: $G=argmin(Fitness(p_{i})$;

12: end if

13: for j:1 to D do do

14: $r_{1}=rand(0,1);r_{2}=rand(0,1)$;

15: $\phi=c_{1}.r_{1}/(c_{1}r_{1}+c_{2}r_{2})$;

16: $p_{i,j}=\phi.P_{i,j}+(1-\phi).G_{j}$;

17: if $rand(0,1)<0.5$ then

18: $x_{i,j}=p_{i,j}+\alpha.abs(mbest_{i,j}-x_{i,j}).log(1/u)$

19: else

20: $x_{i,j}=p_{i,j}-\alpha.abs(mbest_{i,j}-x_{i,j}).log(1/u)$

21: end if

22: end for

23: $Fx(i)$=Fitness$(x(i,:))$;

24: end for

25: $[SF,Sx]=sort(Fx)$;

26: $Z=round((N-1)/S)$; $\triangleright$ S is selection parameter

27: $x(Sx((N-Z+1):N))=x(Sx(1:Z))$; $\triangleright$ Sort the particles from
best to worst position

28:end for

29:Until termination criterion is met

The natural selection method is used to enhance the convergence
characteristics of proposed QPSO-CD algorithm, where the fitter solutions are
used for the next iteration. The procedure of selection method for N particles
is as follows:

$F(X(t))=\\{F(x_{1}(t)),F(x_{2}(t)),...,F(x_{N}(t))\\}$ (6)

where $X(t)$ is position vector of particles at time t and $F(X(t))$ is the
fitness function of swarm. Next step is to sort the particles according to
their fitness values from best one to worst position such that

$\begin{split}F(X^{{}^{\prime}}(t))=\\{F(x_{1}^{{}^{\prime}}(t)),F(x_{2}{{}^{\prime}}(t)),...,&F(x_{N}{{}^{\prime}}(t))\\}\\\
X^{{}^{\prime}}(t)=\\{x_{1}^{{}^{\prime}}(t),x_{2}{{}^{\prime}}(t),...,x_{N}{{}^{\prime}}(t)\\}\end{split}$
(7)

In Algorithm 1, $SF$ and $Sx$ are the sorting functions of fitness and
position respectively. On the basis of natural selection parameters and
fitness values, the positions of swarm particles are updated for the next
iteration,

$\begin{split}X^{{}^{\prime}}(t)=\\{x_{1}^{{}^{\prime\prime}}(t),x_{2}{{}^{\prime\prime}}(t),...,x_{S}{{}^{\prime\prime}}(t)\\},&\\\
X^{{}^{\prime\prime}}_{k}(t)=\\{x_{1}^{{}^{\prime}}(t),x_{2}{{}^{\prime}}(t),...,x_{Z}{{}^{\prime}}(t)\\}\end{split}$
(8)

where $(1\leq k\leq S)$, S denotes the selection parameter, Z signifies the
number of best positions selected according to fitness values such that
$S=N/Z$ and $X^{{}^{\prime\prime}}(t)$ is updated position vector of
particles. The selection parameter S is generally set as 2 to replace the half
of worst positions with the half of best positions of particles. It improves
the precision of the direction of particles, protect the global searching
capability and speed up the convergence.

##  V Experimental results

The performance of proposed QPSO-CD algorithm is investigated on
representative benchmark functions, given in Table 1. Further, the results are
compared with classical PSO (PSO), standard QPSO, QPSO with delta potential
(QDPSO) and QPSO with mutation operator (QPSO-MO). The details of numerical
benchmark functions are given in Table 1.

Table 1: Details of benchmark functions Test function | Initial range
---|---
Sphere function
$f_{1}(x)=\sum\limits_{i=1}^{n}x_{i}^{2}$ | (-100, 100)
Rosenbrock function
$f_{2}(x)=\sum\limits_{i=1}^{n}100(x_{i+1}-x_{i}^{2})^{2}+(x_{i}-1)^{2}$ | (-5.12, 5.12)
Greiwank function
$f_{3}(x)=\dfrac{1}{4000}\sum\limits_{i=1}^{n}x_{i}^{2}-\prod_{i=1}^{n}cos(\dfrac{x_{i}}{\sqrt{i+1}})+1$ | (-600, 600)
Rastrigrin function
$f_{4}(x)=\sum\limits_{i=1}^{n}(x_{i}^{2}-10cos(2\pi x_{i})+10)$ | (-5.12, 5.12)

The performance of QPSO has been widely tested for various test functions.
Initially, we have considered four representative benchmark functions to
determine the reliability of QPSO-CD algorithm. For all the experiments, the
size of population is 20, 40 and 80 and dimension sizes are 10, 20 and 30. The
parameters for QPSO-CD algorithm are as follows: the value of $\alpha$
decreases from 1.0 to 0.5 linearly, the natural selection parameter S=2 is
taken, $c_{1},c_{2}$ correlation coefficients are set equal to 2.

The mean best fitness values of PSO, QPSO, QDPSO, QPSO-MO and QPSO-CD are
recorded for 1000, 1500 and 2000 runs of each function. Fig. 2 to Fig. 5
depict the performance of functions $f_{1}$ to $f_{4}$ with respect to mean
best fitness against the number of iterations. In Table 2, P denotes the
population, dimension is represented by D and G stands for generation. The
numerical results of QPSO-CD shown optimal solution with fast convergence
speed and high accuracy. The results shown that QPSO-CD performs better on
Rosenbrock function than its counterparts in some cases. When the size of
population is 20 and dimension is 30, the results of proposed algorithm are
not better than QPSO-MO, but QPSO-CD performs better than PSO, QPSO and QDPSO.
The performance of QPSO-CD is significantly better than its variants on
Greiwank and Rastrigrin functions. It has outperformed other algorithms and
obtained optimal solution (near zero) for Greiwank function. In most of the
cases, QPSO-CD is more efficient and outperformed the other algorithms.

Figure 2: Effectiveness of QPSO-CD for sphere function $f_{1}$

Figure 3: Effectiveness of QPSO-CD for Rosenbrock function $f_{2}$

Figure 4: Effectiveness of QPSO-CD for Greiwank function $f_{3}$

Figure 5: Effectiveness of QPSO-CD for Rastrigrin function $f_{4}$

Table 2: Comparison results of Sphere and Rosenbrock functions | Sphere function | Rosenbrock function
---|---|---
P | D | G | PSO | QPSO | QDPSO | QPSO-MO | QPSO-CD | PSO | QPSO | QDPSO | QPSO-MO | QPSO-CD
20 | 10 | 1000 | 0.0 | 4.01e-40 | 1.513e-49 | 1.508e-48 | 1.738e-50 | 95.10 | 58.41 | 14.22 | 22.18 | 34.67
20 | 1500 | 0.0 | 2.58e-21 | 1.339e-30 | 1.296e-31 | 1.032e-30 | 204.38 | 110.5 | 175.31 | 68.40 | 54.76
30 | 2000 | 0.0 | 2.08e-13 | 1.953e-21 | 1.918e-21 | 1.808e-21 | 314.46 | 148.5 | 242.37 | 113.30 | 122.5
40 | 10 | 1000 | 0.0 | 2.73e-67 | 1.087e-73 | 1.146e-51 | 1.154e-72 | 70.28 | 10.42 | 15.86 | 7.985 | 8.843
20 | 1500 | 0.0 | 4.84e-28 | 1.397e-42 | 1.417e-42 | 1.237e-41 | 178.98 | 48.45 | 112.46 | 52.93 | 41.77
30 | 2000 | 0.0 | 2.02e-25 | 2.850e-30 | 2.471e-28 | 1.946e-23 | 288.58 | 58.32 | 76.42 | 64.19 | 58.04
80 | 10 | 1000 | 0.0 | 7.66e-95 | 5.553e-90 | 4.872e-71 | 6.437e-72 | 36.29 | 8.853 | 36.34 | 5.715 | 7.419
20 | 1500 | 0.0 | 1.62e-60 | 1.654e-54 | 1.677e-58 | 1.609e-62 | 84.78 | 34.88 | 23.54 | 24.45 | 21.78
30 | 2000 | 0.0 | 2.05e-44 | 1.042e-40 | 1.131e-42 | 1.128e-41 | 202.58 | 52.17 | 70.81 | 45.22 | 40.97
Table 3: Comparison results of Greiwank and Rastrigrin functions | Greiwank function | Rastrigrin function
---|---|---
P | D | G | PSO | QPSO | QDPSO | QPSO-MO | QPSO-CD | PSO | QPSO | QDPSO | QPSO-MO | QPSO-CD
20 | 10 | 1000 | 0.089 | 0.078 | 0.1003 | 0.0732 | 0.072 | 5.526 | 5.349 | 4.969 | 4.478 | 4.051
20 | 1500 | 0.0300 | 0.2001 | 0.0086 | 0.0189 | 0.0078 | 23.17 | 21.28 | 17.08 | 15.63 | 13.22
30 | 2000 | 0.0181 | 0.0122 | 0.0544 | 0.0103 | 0.0026 | 46.29 | 32.57 | 48.61 | 27.80 | 31.48
40 | 10 | 1000 | 0.0826 | 0.055 | 0.048 | 0.0520 | 0.041 | 3.865 | 3.673 | 2.032 | 3.383 | 2.100
20 | 1500 | 0.0272 | 0.0149 | 0.0004 | 0.0247 | 0.0106 | 15.68 | 14.37 | 10.94 | 11.01 | 10.77
30 | 2000 | 0.0125 | 0.0117 | 0.0009 | 0.0105 | 0.0102 | 37.13 | 23.01 | 21.37 | 21.01 | 21.19
80 | 10 | 1000 | 0.0723 | 0.0341 | 0.0 | 0.0542 | 0.0702 | 2.562 | 2.234 | 0.923 | 2.183 | 1.943
20 | 1500 | 0.0274 | 0.0189 | 0.0 | 0.0194 | 0.0161 | 12.35 | 9.66 | 6.955 | 8.075 | 7.021
30 | 2000 | 0.0123 | 0.0118 | 0.0 | 0.0082 | 0.0031 | 26.89 | 17.48 | 18.13 | 14.99 | 11.73

##  VI Correctness and Time Complexity Analysis of a QPSO-CD Algorithm

In this Section, the correctness and time complexity of a proposed algorithm
QPSO-CD is analyzed and compared with the classical PSO algorithm.

######  Theorem 1.

The sequence of random variables $\\{S_{n},n\geq 0\\}$ generated by QPSO with
Cauchy distribution converges to zero in probability as n approaches infinity.

###### Proof.

Recall, the probability density function of standard Cauchy distribution and
its convergence probability rudolph1997local  are given as

$\begin{gathered}f(s)=\dfrac{1}{\pi(1+s^{2})}~{}\text{for}-\infty<s<\infty,\\\
P(x\leq S_{n}\leq
y)=\dfrac{1}{\pi}\int_{y}^{x}\dfrac{ds}{(1+s^{2})},~{}\forall x\leq
y\end{gathered}$ (9)

Consider a random variable $Q_{n}$ interpreted as

$Q_{n}=\alpha S_{n},~{}\alpha=\dfrac{1}{n^{\lambda}}$

where $\lambda$ denotes a fixed positive constant. Correspondingly, the
probability density function can be calculated as

$\displaystyle P(Q_{n}\leq q)$ $\displaystyle=P(\alpha S_{n}\leq q)$
$\displaystyle=P\left(S_{n}\leq\dfrac{q}{\alpha}\right)$
$\displaystyle=\int_{-\infty}^{\dfrac{q}{\alpha}}f(s)ds.\dfrac{dP(q_{n}\leq
q)}{dq}$ $\displaystyle=\dfrac{1}{\alpha}f\left(\dfrac{q}{\alpha}\right)$

i.e. the probability density function of random variable $Q_{n}$.

$\displaystyle P(|Q_{n}|>\xi)$ $\displaystyle=P(|S_{n}|>\xi n^{\lambda})$
$\displaystyle=P(S_{n}>\xi n^{\lambda})+P(S_{n}>-\xi n^{\lambda})$
$\displaystyle=P(\xi n^{\lambda}<S_{n}<\infty)+P(-\infty<S_{n}<-\xi
n^{\lambda})$

Using Eq. (9), the probability density function of random variable $Q_{n}$
becomes

$\displaystyle P(|Q_{n}|>\xi)$ $\displaystyle=\dfrac{1}{\pi}\int_{\xi
n^{\lambda}}^{\infty}\dfrac{ds}{\pi(1+s^{2})}+\dfrac{1}{\pi}\int_{\infty}^{-\xi
n^{\lambda}}\dfrac{ds}{\pi(1+s^{2})}$
$\displaystyle=\left[1+\dfrac{1}{\pi}\int_{\xi n^{\lambda}}^{-\xi
n^{\lambda}}\dfrac{ds}{\pi(1+s^{2})}\right]$
$\displaystyle=1-\dfrac{1}{\pi}\int_{-\xi n^{\lambda}}^{\xi
n^{\lambda}}\dfrac{ds}{\pi(1+s^{2})}=0~{}\text{as}~{}n\rightarrow\infty$

This completes the proof of the theorem. ∎

######  Definition 1.

Let $\\{S_{n}\\}$ a random sequence of variables. It is converges to some
random variable s with probability 1, if for every $\xi>0$ and $\lambda>0$,
there exists $n_{1}(\xi,\lambda)$ such that
$P(|S_{n}-s|<\xi)>1-\lambda,\forall n>n_{1}$ or

$P\left(\lim_{n\rightarrow\infty}|S_{n}-s|<\xi\right)=1$ (10)

The efficiency of the QPSO-CD algorithm is evaluated by number of steps needed
to reach the optimal region $R(\xi)$. The method is to evaluate the
distribution of number of steps needed to hit $R(\xi)$ by comparing the
expected value and moments of distribution. The total number of stages to
reach the optimal region is determined as $W(\xi)=inf\\{n\mid f_{n}\in
R(\xi)\\}$. The variance $V(W(\xi))$ and expectation value $E(W(\xi))$ are
determined as

$E(W(\xi))=\sum_{n=0}^{\infty}nx_{n},$ (11)

$\begin{split}V(W(\xi))&=E(W^{2}(\xi))-\\{E(W(\xi))\\}^{2}\\\
&=\sum_{n=0}^{\infty}n^{2}x_{n}-\left(\sum_{n=0}^{\infty}nx_{n}\right)^{2}\end{split}$
(12)

In fact, the $E(W(\xi))$ depends upon the convergence of
$\sum_{n=0}^{\infty}nx_{n}$. It is needed that $\sum_{j=0}^{\infty}x_{n}=1$,
so that QPSO-CD can converge globally. The number of objective function
evaluations are used to measure time. The main benefit of this approach is
that it shows relationship between processor and measure time as the
complexity of objective function increases. We used Sphere function
$f(x)=x^{T}.x$ with a linear constraint $g(x)=\sum^{n}_{j=0}x_{j}\geq 0$ to
compute the time complexity. It has minimum value at 0. The value of optimal
region is set as $R(\xi)=R(10^{-4}).$ To determine the time complexity, the
algorithms PSO and QPSO-CD are executed 40 times on $f(x)$ with initial scope
[-10, 10]N, where N denotes the dimension. We determine the mean number of
objective function evaluations ($W(\xi)$), the variance ($V(W(\xi))$), the
standard deviation (SD) ($\sigma_{W(\xi)}$), the standard error (SE)
($\sigma_{W(\xi)}/\sqrt{40}$) and ratio of mean and dimension ($W(\xi)/N$).
The contraction coefficient $\alpha=0.75$ is used for QPSO-CD and constriction
coefficient $\chi=0.73$ for PSO with acceleration factors $c_{1}=c_{2}$=2.25.

Table 4: Results of the time complexity for QPSO-CD algorithm Dimension (N) | Mean | Variance | SD | SE | Mean/N
---|---|---|---|---|---
2 | 302.38 | 4164.3 | 64.532 | 10.203 | 151.19
3 | 452.18 | 4541.9 | 67.394 | 10.655 | 150.72
4 | 621.29 | 5208.2 | 72.168 | 11.410 | 155.32
5 | 755.88 | 6675.0 | 81.701 | 12.918 | 151.17
6 | 879.13 | 8523.7 | 92.324 | 14.597 | 146.52
7 | 1022.06 | 9575.4 | 97.854 | 15.472 | 146.00
8 | 1158.52 | 10269.7 | 101.341 | 16.023 | 144.81
9 | 1308.17 | 12053.4 | 109.788 | 17.359 | 145.35
10 | 1459.3 | 12648.3 | 112.465 | 17.782 | 145.93
Table 5: Results of the time complexity for PSO algorithm Dimension (N) | Mean | Variance | SD | SE | Mean/N
---|---|---|---|---|---
2 | 691.4 | 17297.5 | 131.52 | 20.795 | 345.7
3 | 979.1 | 22281.5 | 149.27 | 23.601 | 326.3
4 | 1167.2 | 24282.9 | 155.72 | 24.638 | 291.8
5 | 1328.7 | 21853.7 | 147.83 | 23.373 | 265.7
6 | 1489.9 | 32008.7 | 178.91 | 28.288 | 248.3
7 | 1744.3 | 502297 | 224.12 | 35.436 | 249.1
8 | 1978.5 | 41233.3 | 203.06 | 31.106 | 247.3
9 | 2259.1 | 36217.8 | 190.31 | 30.090 | 251.0
10 | 2604.2 | 43559.8 | 208.71 | 32.999 | 260.4
Figure 6: Time complexity results for PSO and QPSO-CD Figure 7: Comparison of
Correlation coefficient of PSP and QPSO-CD

Table 1 and 2 show the statistical results of time complexity test for QPSO-CD
and PSO algorithm, respectively. Fig 6 indicates that the time complexity of
proposed algorithm increases non-linearly as the dimension increases. However,
the time complexity of PSO algorithm increases adequately linearly. Thus, the
time complexity of QPSO-CD is lower than PSO algorithm. In Fig 7, QPSO-CD
shows a strong correlation between $W(\xi)$ and N, i.e. the correlation
coefficient is 0.9996. For PSO, the linear correlation coefficient is 0.9939,
which is not so phenomenal as that in case of QPSO-CD. The relationship
between mean and dimension clearly shows that the value of correlation
coefficient is fairly stable for QPSO-CD as compared to PSO algorithm.

##  VII QPSO-CD for Constraint Engineering Design Problems

There exists several approaches for handling constrained optimization
problems. The basic principle is to convert the constrained optimization
problem to unconstrained by combining objective function and penalty function
approach. Further, minimize the newly formed objective function with any
unconstrained algorithm. Generally, the constrained optimization problem can
be described as in Eq. (13).

The objective is to minimize the objective function f(x) subjected to equality
$(h_{j}(x))$ and inequality $(g_{i}(x))$ constrained functions, where $p(i)$
is the upper bound and $q(i)$ denotes the search space lower bound. The strict
inequalities of form $g_{i}(x)\geq 0$ can be converted into $-g_{i}(x)\leq 0$
and $h_{i}(x)$ equality constraints can be converted into inequality
constraints $h_{i}(x)\geq 0$ and $h_{i}(x)\leq 0$. Sun et al. sun2007using
adopted non-stationary penalty function to address non-linear programming
problems using QPSO. Coelho dos2010gaussian  used penalty function with some
positive constant i.e. set to 5000. We adopted the same approach and replace
the constant with dynamically allocated penalty value.

$\begin{split}&\min\limits_{x}=f(x)\\\ \text{subject to}\\\ &g_{i}(x)\leq
0,i=0,1,...n-1\\\ &h_{j}(x)=0,j=1,2,...r\\\ &p(i)\leq x_{i}\leq q(i),1\leq
i\leq m\\\ &x=\\{x_{1},x_{2},x_{3},...,x_{m}\\}\end{split}$ (13)

Usually, the procedure is to find the solution for design variables that lie
in search space upper and lower bound constraints such that
$x_{i}\in[p(i),q(i)]$. If solution violates any of the constraint, then the
following rules are applied

$\begin{split}&x_{i}=x_{i}+\\{p(x_{i})-q(x_{i})\\}.~{}rand[0,1]\\\
&x_{i}=x_{i}-\\{p(x_{i})-q(x_{i})\\}.~{}rand[0,1]\end{split}$ (14)

where rand[0, 1] is randomly distributed function to select value between 0
and 1. Finally, the unconstrained optimization problem is solved using
dynamically modified penalty values according to inequality constraints
$g_{i}(x)$. Thus, the objective function is evaluated as

$F(x)=\left\\{\begin{array}[]{ll}f(x)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if}~{}g_{i}(x)\leq
0\\\ f(x)+y(t).\sum\limits_{i=1}^{n}g_{i}(x)~{}~{}\text{if}~{}g_{i}(x)>0\\\
\end{array}\right\\}$ (15)

where $f(x)$ is the main objective function of optimization problem in Eq.
(13), t is the iteration number and $y(t)$ represents the dynamically
allocated penalty value.

In this Section, QPSO-CD is tested for three-bar truss, tension/compression
spring and pressure vessel design problems consisting different members and
constraints. The performance of QPSO-CD is compared and analyzed with the
results of PSO, QPSO, and SP-QPSO algorithms as reported in the literature.

###  VII.1 Three-bar truss design problem

Three-bar truss is a constraint design optimization problem, which has been
widely used to test several methods. It consists cross-section areas of three
bars $x_{1}$ (and $x_{3}$) and $x_{2}$ as design variables. The aim of this
problem is to minimize the weight of truss subject to maximize the stress on
these bars. The structure should be symmetric and subjected to two constant
loadings $P_{1}=P_{2}=P$ as shown in Fig 8. The mathematical formulation of
two design bars ($x_{1}$, $x_{2}$) and three restrictive mathematical
functions are described as:

$\begin{split}\text{Min}~{}f(x)&=(2\sqrt{2}x_{1}+x_{2}).l\\\ \text{Subject
to:}&\\\
g_{1}(x)&=\dfrac{\sqrt{2}x_{1}+x_{2}}{\sqrt{2}x_{1}^{2}+2x_{1}x_{2}}P-\sigma\leq
0,\\\ g_{2}(x)&=\dfrac{x_{2}}{\sqrt{2}x_{1}^{2}+2x_{1}x_{2}}P-\sigma\leq 0\\\
\end{split}$ (16) $g_{3}(x)=\dfrac{1}{x_{1}+\sqrt{2}x_{2}}P-\sigma\leq
0,~{}\text{where}$ $0\leq x_{1},x_{2}\leq 1,l=100cm,P=\sigma=2KN/cm^{2}$

Figure 8: Structure of Three-bar truss Table 6: Comparison of optimal results for three-bar truss problem Variables | PSO | QPSO | SP-QPSO | QPSO-CD
---|---|---|---|---
$x_{1}$ | 0.78911058 | 0.788649 | 0.788796 | 0.788658
$x_{2}$ | 0.40702683 | 0.408322 | 0.407898 | 0.40828488
$g_{1}(x)$ | -6.6720e-06 | 1.6313e-07 | 6.4748e-06 | 9.00037e-06
$g_{2}(x)$ | -1.4655 | -1.4640 | -1.4644 | -1.4640
$g_{3}(x)$ | -0.5345 | -0.5359 | -0.5354 | -0.5359
$f(x)$ | 263.89686 | 263.89584 | 263.89500 | 263.89465

The results are obtained by QPSO-CD are compared with its counterparts in
Table 6. For three-bar truss problem, QPSO-CD is superior to optimal solutions
previously obtained in literature. The difference of best solution obtained by
QPSO-CD among other algorithms is shown in Fig 9.

Figure 9: Optimal results of PSO, QPSO, SP-QPSO and QPSO-CD algorithms for
three bar truss problem

###  VII.2 Tension/Compression spring design problem

The main aim is to lessen the volume V of a spring subjected to tension load
constantly as shown in Fig 10. Using the symmetry of structure, there are
practically three design variables ($x_{1},x_{2},x_{3}$), where $x_{1}$ is the
wire diameter, the coil diameter is represented by $x_{2}$ and $x_{3}$ denotes
the total number of active coils. The mathematical formulation for this
problem is described as:

$\begin{split}\text{Min}~{}f(x)&=(x_{3}+2)x_{2}x_{1}^{2},\\\ \text{Subject
to:}&\\\ g_{1}(x)&=\dfrac{1-x_{2}^{3}x_{3}}{71785x_{1}6{4}}\leq 0,\\\
g_{2}(x)&=\dfrac{4x_{2}^{2}-x_{1}x_{2}}{12566(x_{2}x_{1}^{3}-x_{1}^{4})}+\dfrac{1}{5108x_{1}^{2}}-1\leq
0,\\\ g_{3}(x)&=\dfrac{1-140.45x_{1}}{x_{2}^{2}x_{3}}\leq 0,\\\
g_{4}(x)&=\dfrac{x_{2}+x_{1}}{1.5}-1\leq 0,\\\ \text{where}&\\\ 0.05\leq
x_{1}&\leq 2,0.25\leq x_{2},\leq 1.3,2\leq x_{3}\leq 15,\\\ \end{split}$ (17)

Figure 10: Structure of tension/compression spring Table 7: Comparison of optimal results for tension spring design problem Variables | PSO | QPSO | SP-QPSO | QPSO-CD
---|---|---|---|---
$x_{1}$ | 0.0516 | 0.0524 | 0.05 | 0.0513
$x_{2}$ | 0.3542 | 0.2505 | 0.25 | 0.2502
$x_{3}$ | 11.7942 | 2 | 2 | 2
$g_{1}(x)$ | -2.3006e-02 | 0.93145095 | 0.93034756 | 4.11004e-06
$g_{2}(x)$ | -5.6059e-03 | -0.17471558 | -0.16568318 | -0.17352479
$g_{3}(x)$ | -3.9057 | -50.67 | -48.180 | -49.561
$g_{4}(x)$ | -0.7294 | -0.79986567 | -0.80 | -0.799
$f(x)$ | 0.01305 | 0.00275 | 0.00250 | 0.00263

It has been observed that QPSO algorithm with Cauchy distribution and natural
selection strategy is robust and obtains optimal solutions than PSO and QPSO,
shown in Table 7. The difference between best solutions found by QPSO-CD
($f(x)=0.00263$) and other algorithms for tension spring design problem are
reported in Fig 11.

Figure 11: Results of PSO, QPSO, SP-QPSO and QPSO-CD methods for tension
spring design problem

###  VII.3 Pressure vessel design problem

Initially, Kannan and Kramer kannan1994augmented  studied the pressure vessel
design problem with the main aim to reduce the total fabricating cost.
Pressure vessels can be of any shape. For engineering purposes, a cylindrical
design capped by hemispherical heads at both ends is widely used
sandgren1990nonlinear . Fig 12 describes the structure of pressure vessel
design problem. It consists four design variables ($x_{1},x_{2},x_{3},x_{4}$),
where $x_{1}$ denotes the shell thickness $(T_{s})$, $x_{2}$ is used for head
thickness ($T_{h}$), $x_{3}$ denotes the inner radius (R) and $x_{4}$
represents the length of vessel (L). The objective function and constraint
equations are described as:

Figure 12: Design of pressure vessel

$\begin{split}\text{Min}~{}f(x)&=0.6224x_{1}x_{3}x_{4}+1.7781x_{2}x_{3}^{2}+3.166x_{1}^{2}x_{4}\\\
&+19.84x_{1}^{2}x_{3},\\\ \text{Subject to:}&\\\
g_{1}(x)&=-x_{1}+0.0193x_{3}\leq 0,\\\ g_{2}(x)&=-x_{2}+0.00954x_{3}\leq 0,\\\
g_{3}(x)&=-\pi x_{3}^{2}x_{4}-\dfrac{4}{3}\pi x_{3}^{3}+1296000\leq 0,\\\
g_{4}(x)&=x_{4}-240\leq 0,\\\ \text{where}&\\\ 1\times 0.0625&\leq
x_{1},x_{2}\leq 99\times 0.0625,10\leq x_{3},x_{4}\leq 200\\\ \end{split}$
(18)

Table 8: Comparison of optimal results for Pressure vessel design problem Variables | PSO | QPSO | SP-QPSO | QPSO-CD
---|---|---|---|---
$x_{1}$ | 0.8125 | 0.7783 | 0.7782 | 0.7776
$x_{2}$ | 0.4375 | 0.3849 | 0.3845 | 0.3848
$x_{3}$ | 42.0984 | 40.3289 | 40.3206 | 40.3278
$x_{4}$ | 176.6365 | 199.8899 | 199.9988 | 199.8865
$g_{1}(x)$ | -4.500e-15 | 4.777e-05 | -1.242e-05 | 7.2654e-04
$g_{2}(x)$ | -0.035880 | -1.62294e-04 | 1.58523e-04 | -7.2787e-05
$g_{3}(x)$ | -1.164e-10 | -97.39720071 | -63.63686942 | -0.734359
$g_{4}(x)$ | -63.3634 | -40.1100 | -40.0012 | -40.1135
$f(x)$ | 6059.714 | 5886.189 | 5885.268 | 5886.137
Figure 13: Optimal results of PSO, QPSO, SP-QPSO and QPSO-CD techniques for
pressure vessel design problem

The optimal results of QPSO-CD is compared with the SP-QPSO, QPSO and PSO best
results noted in the previous work, and are given in Table 8. The best
solution obtained from QPSO-CD is better than other algorithms as shown in Fig
13.

##  VIII Conclusion

In this paper, a new hybrid quantum particle swarm optimization algorithm is
proposed with natural selection method and Cauchy distribution. The
performance of the proposed algorithm is experimented on four benchmark
functions and the optimal results are compared with existing algorithms.
Further, the QPSO-CD is applied to solve engineering design problems. The
efficiency of QPSO-CD is successfully presented with superiority than
preceding results for three engineering design problems: three-bar truss,
tension/compression spring and pressure vessel. The efficiency of QPSO-CD
algorithm is evaluated by number of steps needed to reach the optimal region
and proved that time complexity of proposed algorithm is lower in comparison
to classical PSO. In the context of convergence, the experimental outcomes
shown that the QPSO-CD converge to get results closer to the superior
solution.

## Additional information

Competing interests: The authors declare no competing interests.

## Acknowledgement

S.Z. acknowledges support in part from the National Natural Science Foundation
of China (Nos. 61602532), the Natural Science Foundation of Guangdong Province
of China (No. 2017A030313378), and the Science and Technology Program of
Guangzhou City of China (No. 201707010194).

## References

  * (1) N. Bohr, et al., The quantum postulate and the recent development of atomic theory, Vol. 3, Printed in Great Britain by R. & R. Clarke, Limited, 1928.
  * (2) H. P. Robertson, The uncertainty principle, Physical Review 34 (1) (1929) 163.
  * (3) L. Wessels, Schrödinger’s route to wave mechanics, Studies in History and Philosophy of Science Part A 10 (4) (1979) 311–340.
  * (4) R. P. Feynman, Simulating physics with computers, International journal of theoretical physics 21 (6-7) (1982) 467–488.
  * (5) J. Wang, Handbook of Finite State Based Models and Applications, CRC press, 2012\.
  * (6) M. Dorigo, G. Di Caro, Ant colony optimization: a new meta-heuristic, in: Proceedings of the 1999 congress on evolutionary computation-CEC99 (Cat. No. 99TH8406), Vol. 2, IEEE, 1999, pp. 1470–1477.
  * (7) A. Colorni, M. Dorigo, V. Maniezzo, et al., Distributed optimization by ant colonies, in: Proceedings of the first European conference on artificial life, Vol. 142, Cambridge, MA, 1992, pp. 134–142.
  * (8) J. Kennedy, R. Eberhart, Particle swarm optimization (pso), in: Proc. IEEE International Conference on Neural Networks, Perth, Australia, 1995, pp. 1942–1948.
  * (9) S.-C. Wang, M.-F. Yeh, A modified particle swarm optimization for aggregate production planning, Expert Systems with Applications 41 (6) (2014) 3069–3077.
  * (10) M. R. AlRashidi, M. E. El-Hawary, A survey of particle swarm optimization applications in electric power systems, IEEE transactions on evolutionary computation 13 (4) (2009) 913–918.
  * (11) A. R. Yıldız, A novel particle swarm optimization approach for product design and manufacturing, The International Journal of Advanced Manufacturing Technology 40 (5-6) (2009) 617.
  * (12) N. A. Latiff, C. C. Tsimenidis, B. S. Sharif, Energy-aware clustering for wireless sensor networks using particle swarm optimization, in: 2007 IEEE 18th International Symposium on Personal, Indoor and Mobile Radio Communications, IEEE, 2007, pp. 1–5.
  * (13) S.-W. Lin, K.-C. Ying, S.-C. Chen, Z.-J. Lee, Particle swarm optimization for parameter determination and feature selection of support vector machines, Expert systems with applications 35 (4) (2008) 1817–1824.
  * (14) Y. Wang, S. D. Mohanty, Particle swarm optimization and gravitational wave data analysis: Performance on a binary inspiral testbed, Physical Review D 81 (6) (2010) 063002.
  * (15) M. E. Normandin, S. D. Mohanty, T. S. Weerathunga, Particle swarm optimization based search for gravitational waves from compact binary coalescences: Performance improvements, Physical Review D 98 (4) (2018) 044029.
  * (16) J. Sun, B. Feng, W. Xu, Particle swarm optimization with particles having quantum behavior, in: Proceedings of the 2004 congress on evolutionary computation (IEEE Cat. No. 04TH8753), Vol. 1, IEEE, 2004, pp. 325–331.
  * (17) A. Narayanan, M. Moore, Quantum-inspired genetic algorithms, in: Proceedings of IEEE international conference on evolutionary computation, IEEE, 1996, pp. 61–66.
  * (18) M. Yuanyuan, L. Xiyu, Quantum inspired evolutionary algorithm for community detection in complex networks, Physics Letters A 382 (34) (2018) 2305–2312.
  * (19) J. Liu, W. Xu, J. Sun, Quantum-behaved particle swarm optimization with mutation operator, in: 17th IEEE international conference on tools with artificial intelligence (ICTAI’05), IEEE, 2005, pp. 4–pp.
  * (20) V. Protopopescu, J. Barhen, Solving a class of continuous global optimization problems using quantum algorithms, Physics Letters A 296 (1) (2002) 9–14.
  * (21) A. S. Bhatia, M. K. Saggi, A. Kumar, S. Jain, Matrix product state–based quantum classifier, Neural computation 31 (7) (2019) 1499–1517.
  * (22) A. S. Bhatia, M. K. Saggi, Implementing entangled states on a quantum computer, arXiv preprint: 1811.09833 (2018).
  * (23) J. Sun, W. Xu, B. Feng, A global search strategy of quantum-behaved particle swarm optimization, in: IEEE Conference on Cybernetics and Intelligent Systems, 2004., Vol. 1, IEEE, 2004, pp. 111–116.
  * (24) J. Sun, J. Liu, W. Xu, Using quantum-behaved particle swarm optimization algorithm to solve non-linear programming problems, International Journal of Computer Mathematics 84 (2) (2007) 261–272.
  * (25) J. Sun, W. Fang, V. Palade, X. Wu, W. Xu, Quantum-behaved particle swarm optimization with gaussian distributed local attractor point, Applied Mathematics and Computation 218 (7) (2011) 3763–3775.
  * (26) L. dos Santos Coelho, Gaussian quantum-behaved particle swarm optimization approaches for constrained engineering design problems, Expert Systems with Applications 37 (2) (2010) 1676–1683.
  * (27) Y. Li, R. Xiang, L. Jiao, R. Liu, An improved cooperative quantum-behaved particle swarm optimization, Soft Computing 16 (6) (2012) 1061–1069.
  * (28) Y. Peng, Y. Xiang, Y. Zhong, Quantum-behaved particle swarm optimization algorithm with lévy mutated global best position, in: 2013 Fourth International Conference on Intelligent Control and Information Processing (ICICIP), IEEE, 2013, pp. 529–534.
  * (29) H. Ali, W. Shahzad, F. A. Khan, Energy-efficient clustering in mobile ad-hoc networks using multi-objective particle swarm optimization, Applied Soft Computing 12 (7) (2012) 1913–1928.
  * (30) A. S. Bhatia, R. K. Cheema, Analysing and implementing the mobility over manets using random way point model, International Journal of Computer Applications 68 (17) (2013) 32–36.
  * (31) Z. Zhisheng, Quantum-behaved particle swarm optimization algorithm for economic load dispatch of power system, Expert Systems with Applications 37 (2) (2010) 1800–1803.
  * (32) J. Sun, J. Liu, W. Xu, Qpso-based qos multicast routing algorithm, in: Asia-Pacific Conference on Simulated Evolution and Learning, Springer, 2006, pp. 261–268.
  * (33) M. Geis, M. Middendorf, Particle swarm optimization for finding RNA secondary structures, International Journal of Intelligent Computing and Cybernetics (2011).
  * (34) A. S. Bhatia, A. Kumar, Modeling of RNA secondary structures using two-way quantum finite automata, Chaos, Solitons & Fractals 116 (2018) 332–339.
  * (35) A. Bagheri, H. M. Peyhani, M. Akbari, Financial forecasting using anfis networks with quantum-behaved particle swarm optimization, Expert Systems with Applications 41 (14) (2014) 6235–6250.
  * (36) E. Davoodi, M. T. Hagh, S. G. Zadeh, A hybrid improved quantum-behaved particle swarm optimization–simplex method (iqpsos) to solve power system load flow problems, Applied Soft Computing 21 (2014) 171–179.
  * (37) S. Omkar, R. Khandelwal, T. Ananth, G. N. Naik, S. Gopalakrishnan, Quantum behaved particle swarm optimization (qpso) for multi-objective design optimization of composite structures, Expert Systems with Applications 36 (8) (2009) 11312–11322.
  * (38) D. Fatemeh, C. Loo, G. Kanagaraj, Shuffled complex evolution based quantum particle swarm optimization algorithm for mechanical design optimization problems, Journal of Modern Manufacturing Systems and Technology 2 (1) (2019) 23–32.
  * (39) S. Prithi, S. Sumathi, Ld2fa-pso: A novel learning dynamic deterministic finite automata with pso algorithm for secured energy efficient routing in wireless sensor network, Ad Hoc Networks 97 (2020) 102024.
  * (40) A. Ambainis, R. Freivalds, 1-way quantum finite automata: strengths, weaknesses and generalizations, in: Proceedings 39th Annual Symposium on Foundations of Computer Science (Cat. No. 98CB36280), IEEE, 1998, pp. 332–341.
  * (41) A. S. Bhatia, A. Kumar, Quantum finite automata: survey, status and research directions, arXiv preprint:1901.07992 (2019).
  * (42) A. S. Bhatia, A. Kumar, On the power of two-way multihead quantum finite automata, RAIRO-Theoretical Informatics and Applications 53 (1-2) (2019) 19–35.
  * (43) J. Kennedy, Particle swarm optimization, Encyclopedia of machine learning (2010) 760–766.
  * (44) Y. Shi, et al., Particle swarm optimization: developments, applications and resources, in: Proceedings of the 2001 congress on evolutionary computation (IEEE Cat. No. 01TH8546), Vol. 1, IEEE, 2001, pp. 81–86.
  * (45) G. Rudolph, Local convergence rates of simple evolutionary algorithms with cauchy mutations, IEEE Transactions on Evolutionary Computation 1 (4) (1997) 249–258.
  * (46) B. Kannan, S. N. Kramer, An augmented lagrange multiplier based method for mixed integer discrete continuous optimization and its applications to mechanical design, Journal of mechanical design 116 (2) (1994) 405–411.
  * (47) E. Sandgren, Nonlinear integer and discrete programming in mechanical design optimization, Journal of Mechanical Design 112 (2) (1990) 223–229.


# Viscoelastic coarsening of quasi-2D foam

Chiara Guidolin  Université Paris-Saclay, CNRS, Laboratoire de Physique des
Solides, Orsay 91405, France  Jonatan Mac Intyre  Aalto University, Department
of Applied Physics, Espoo 02150, Finland  Emmanuelle Rio  Université Paris-
Saclay, CNRS, Laboratoire de Physique des Solides, Orsay 91405, France  Antti
Puisto  Aalto University, Department of Applied Physics, Espoo 02150, Finland
Anniina Salonen  Université Paris-Saclay, CNRS, Laboratoire de Physique des
Solides, Orsay 91405, France

###### Abstract

Foams are unstable jammed materials. They evolve over timescales comparable to
their "time of use", which makes the study of their destabilisation mechanisms
crucial for applications. In practice, many foams are made from viscoelastic
fluids, which are observed to prolong their lifetimes. Despite their
importance we lack understanding of the coarsening mechanism in such systems.
We probe the effect of continuous phase viscoelasticity on foam coarsening
with foamed emulsions. We show that bubble size evolution is strongly slowed
down and foam structure hugely impacted. The main mechanisms responsible are
the absence of continuous phase redistribution and a non-trivial link between
foam structure and mechanical properties. These combine to give spatially
heterogeneous coarsening. Beyond their importance in the design of foamy
materials, the results give a macroscopic vision of phase separation in a
viscoelastic medium.

Foams coarsen as gas diffuses between bubbles due to differences in Laplace
pressure. The rate of coarsening is set by foam topology through the
organisation of the neighbouring bubbles[1, 2]. The growth law is
theoretically predicted and experimentally verified in two limiting cases of
very dry or very wet foams. In both cases the foam reaches a self-similar
scaling state with invariant statistical distributions [3], which leads to a
power-law growth of the mean bubble size $R$. In dry foams gas diffusion
occurs exclusively through the thin films separating the bubbles and
$R(t)\propto t^{1/2}$, while in dilute bubble dispersions the gas transfer
through the bulk phase is captured within a mean-field description and
$R(t)\propto t^{1/3}$[4].

If gas diffusion is no longer the rate-limiting step in the coarsening
process, the evolution could be different. In wet foams made from viscoelastic
fluids, lamellar[5] or smectic[6], exponents below $1/3$ have been measured.

Bubbles in a viscoelastic continuous phase stop coarsening if the surrounding
material is sufficiently stiff to oppose the growth (or shrinking) of the
bubbles. The conditions for arrest have been predicted for both single bubbles
and foams [7, 8]. Bey et al. worked with monodisperse bubbles to show that,
indeed, if the elastic energy of the medium is higher than the surface energy
coarsening is halted [9]. Lesov et al. measured an arrest once the continuous
phase yield stress was higher than Laplace pressure [10].

The relevant parameter to compare surface tension effects to bulk elasticity
is the elastocapillary number $Ca_{\text{el}}$ [11]. In the case of capillary
inclusions (liquid drops or gas bubbles) in an elastic matrix,
$Ca_{\text{el}}$ is defined as the ratio of the bulk elastic modulus $G_{0}$
to Laplace pressure, thus $Ca_{\text{el}}=\frac{G_{0}}{\gamma/R}$, where
$\gamma$ is the surface tension and $R$ the bubble (drop) size. Depending on
$Ca_{\text{el}}$, the inclusions can either soften or stiffen the material
[12, 13]. In a foam the bubbles cannot be considered as spherical inclusions,
as they are deformed with thin films between them. However $Ca_{\text{el}}$
has been shown to be the relevant parameter to describe foam elasticity due to
the additivity of the bubble and continuous phase contributions [14, 15].

In practice, we encounter foamy materials in a broad range of
$Ca_{\text{el}}$. Between the limits of capillary-controlled foams (e.g. soap
froths) and solid foams, lies a number of materials with intermediate
$Ca_{\text{el}}$. These are foams made from soft solids (creams, pastes, gels)
or precursors of solid foams. In these systems, $Ca_{\text{el}}$ evolves in
time, due to an increase of the bubble size (coarsening, coalescence) or
solidification of the continuous phase, both of which lead to an increase of
$Ca_{\text{el}}$. Understanding bubble size evolution in such systems is
crucial as it is a key control parameter of the material properties [16]. We
thus explore the evolution of the bubble size and foam structure in a foam
made from a model viscoelastic fluid (an emulsion), in which we vary the
initial elastocapillary number from 1 to 8. In this range we expect the
continuous phase to impact the coarsening process without arresting it.

In concentrated emulsions the elastic and viscous moduli can be easily varied
through changing the volume fraction $\phi$ of the dispersed phase[17]. This
is why we use oil-in-water emulsions, with an average drop size of 2 $\mu$m
and $\phi$ between 0.65 to 0.85, which allows us to vary the continuous phase
elastic modulus from 31 to 506 Pa (Fig. 1a). All the studied emulsions are
predominantly elastic, the loss modulus always being less than 10% of the
storage modulus. We use sodium dodecyl sulphate (SDS) at 30 g/L as the
surfactant to stabilise both the oil/water and the gas/water interfaces. The
choice of such a high SDS concentration ensures that all the interfaces are
well covered.

We choose to work with quasi-2D foams, which consist of a single layer of
bubbles sandwiched between two glass plates (Fig. 1d). Such systems have been
widely used by the foam community to explore many aspects of foam
stability[18, 19] and rheology[20]. The coarsening laws are the same as in 3D
foams, but the specific configuration allows the tracking of individual
bubbles for a fine analysis of the foam structure[21]. Instead of creating the
foam directly between the plates, we first make a 3D foamed emulsion with
small polydisperse bubbles and a liquid (emulsion) fraction $\epsilon\approx
11\%$ using a mixer to incorporate air in the emulsion (Fig. 1c). Then we fill
and close the cell to leave the bubbles to grow. This allows us to start with
quasi-2D foams with small, polydisperse bubbles (Fig. 1d). The images are
skeletonised to measure the bubble size (Fig. 1e). We use $\langle R\rangle$ =
0.5 mm to calculate the initial $Ca_{\text{el}}$, which for the different
samples increases from 0.5 to 8 with increasing $\phi$ (Fig. 1a).

Figure 1: Foamed emulsion setup. a, Storage (down triangles) and loss moduli
(up triangles) of the emulsions used as the foam continuous phases as a
function of oil volume fraction $\phi$. The dashed line shows the expected
scaling of $G^{\prime}\approx\phi(\phi-\phi^{*})$ with $\phi^{*}=0.635$ from
[17], which describes the data well. In the following we use $G_{0}$ for
$G^{\prime}$ for the emulsion elastic modulus. The squares are
$Ca_{\text{el}}$ calculated using $\langle R\rangle$ = 0.5 mm and a surface
tension of 30 mN/m and the solid line guides to $Ca_{\text{el}}=1$. b,
Schematic drawing of the sample generation from 3D foam to the quasi-2D foam.
The gap between the plates is 1 mm. c, Photograph of a foam after generation
seen under the microscope. The edge size is 830 $\mu$m. d, Perspective view of
a quasi-2D foam. e, Top view of a foam with its skeleton overlaid in red. The
edge size is 30 mm.

Photographs of three foams at different oil fractions, taken at
$t=t_{\text{2D}}$, are shown in Fig. 2a. The pictures at $\phi=0.65$ and
$0.75$ are similar and the typical cellular foam structure can be recognised.
But at $\phi=0.85$ the foam already has some peculiar features, such as
elongated bubble shapes, which we will return to later on. Despite the
morphological differences, the bubble size distributions at this point are
very similar in all the foams considered. This can be seen from the
distribution functions of the dimensionless radius $R/\langle R\rangle$ (Fig.
2b), which are also comparable to the steady state distributions measured in
aqueous quasi-2D foams [22].

The temporal evolution of the average radius normalised by the average radius
at $t=t_{\text{2D}}$ of the quasi-2D foams changes considerably with $\phi$
(Fig. 2c). The foam made from the least elastic emulsion at $\phi=0.65$
evolves the most rapidly and with an apparent power law. However, the exponent
is smaller than the expected $t^{1/2}$ for dry foams and closer to $t^{1/3}$
expected for very wet foams, although the foam is rather dry. This is a
consequence of finite size effects observed in moderately dry quasi-2D foams,
where the vertical film size decreases as the foam coarsens at constant liquid
fraction [23]. Even though our systems are dry enough for adjacent bubbles to
share thin films, all samples have coarsening rates equal to or below the
prediction for dilute bubble dispersions. This is the lowest possible rate for
diffusion-limited growth, and any change from diminished gas transport

in the emulsion would only shift the curves, but not change the power laws
observed. The foams made from emulsions at higher $\phi$ evolve very slowly,
and with continuously slowing coarsening rates with increasing bubble size. As
the liquid fraction in the foams is very similar, the differences between the
samples are due to the emulsion.

Figure 2: Initial foams and size evolution. a, Photographs of foams made from
the different emulsions at $t=t_{\text{2D}}$. From top to bottom: $\phi$=0.85,
0.75, 0.65. The edge of each frame is 40 mm. b, Bubble size distributions at
$t=t_{\text{2D}}$ for the five different foams studied, on the left axis the
probability distribution function and on the right axis the cumulative
distribution function. The initial distributions of the foams are almost
identical, and the cumulative distributions overlap almost perfectly. c,
Normalised mean bubble size growth. The black solid line is a $1/3$ powerlaw,
which is the lower limit for diffusion limited coarsening. Inset: Temporal
evolution of the total number of measured bubbles, where the limit for
diffusive ripening is given by the $-2/3$ law.

We have a high number of bubbles $N$ throughout the coarsening process (inset
of Fig. 2c), which allows us to measure bubble size distributions in time with
good statistics. We thus compare the evolution of the dimensionless bubble
size distributions in the foams made from the different emulsions (Fig. 3a).
It is clear that the elasticity of the continuous phase has an impact on the
distributions, which unlike aqueous foams, do not head towards self-similar
distributions obtained at steady state[22]. Although the exact evolution
depends on $\phi$, in all of them a shift of the peak to smaller $R/\langle
R\rangle$ can be observed This shift implies an accumulation of smaller
bubbles along time, which indicates a delay in their disappearance that is at
the origin of the slow coarsening evolutions observed.

We choose to quantify the change in the shape of the distributions through its
third central moment $\mu_{3}^{R}$, which is sensitive to the distribution
asymmetry or skew and, in this case, to the presence of small bubbles. By
plotting $\mu_{3}^{R}$ as a function of the average bubble radius (Fig. 3b),
we can see that $\mu_{3}^{R}$ is initially around 0.03 for all the foams, but
as the bubble size increases it departs from this plateau. The higher $\phi$,
the smaller the bubble size at which the skew deviates. We note that at $\phi$
= 0.65 a deviation is observed as the bubbles reach around 2 mm, but this is
due to the onset of coalescence and hence the physical origin is different
than in the other samples. We define the radius at which $\mu_{3}^{R}$ has
increased by 50% from its average value (from 0.036 to 0.053) as
$R_{\text{sk}}$. This gives us a characteristic radius at which the foam
evolution has deviated from that of classical foams. $R_{\text{sk}}$ decreases
with increased $\phi$, and hence higher $G_{0}$ (Fig. 3c).

Figure 3:  Evolution of bubble size distributions. a, Time evolution of the
dimensionless bubble size distributions in the foams made from emulsions at
different $\phi$. The curves are plotted with a grey scale proportional to the
foam age, where black corresponds to the end of the image acquisition $t_{f}$.
b, Third moment $\mu_{3}^{R}$ versus the mean bubble radius. The solid black
line indicates the average value at early times. The grey zones show the $\pm
50\%$ and $+100\%$ values, which are used to estimate $R_{\text{sk}}$ in d.
The $\mu_{3}^{R}$ have black stars, which indicate the values at $t=96$ hours.
c, Impact of emulsion elastic modulus on $R_{\text{sk}}$. The black line shows
the prediction with $Ca_{\text{el}}^{\text{FE}}=0.54$ and the grey region the
impact of changing liquid fraction on it. d, Evolution of the expected foam
elastic modulus $G_{\text{FE}}$ (solid lines) and of the foam elastocapillary
number $Ca_{\text{el}}^{\text{FE}}$ (dashed lines) as the mean bubble size
grows because of coarsening, for the different $\phi$.

We expect the departure from standard foam coarsening to happen once the
elasticity of the continuous phase becomes large compared to Laplace pressure.
However, in a foam the bubbles are not surrounded by a continuous emulsion
phase of constant elastic modulus, but by other bubbles forming a foam. The
elastic modulus of a 3D foam with an elastic continuous phase can be described
by an additive model combining bubble $G_{F}$ and continuous phase $G_{E}$
contributions as $G_{\text{FE}}=G_{F}+G_{E}$[14]. The foam elasticity is
modeled by $G_{F}=1.6\dfrac{\gamma}{\langle
R\rangle}(1-\epsilon)(0.36-\epsilon)\equiv f(\epsilon)\frac{\gamma}{\langle
R\rangle}$ [24], and the continuous phase by $G_{E}=\epsilon^{2}G_{0}$ [16],
where $\epsilon$ is the foam liquid fraction. Despite the structural
differences between quasi-2D and 3D foams we use the 3D foam description to
capture the dependency of $G_{FE}$ on $R$ and $G_{0}.$ For simplicity we do
not use the coupling function introduced by Gorlier et al[14], as it only has
a weak effect on our estimated moduli and as the function could be system
dependent. The elastic moduli of foams $G_{\text{FE}}$ calculated with the
foam liquid fraction $\epsilon$ and the elastic moduli $G_{0}$ corresponding
to our emulsions, are shown as a function of the average bubble radius in Fig.
3d. At small $R$ the $G_{\text{FE}}$ are all close as they are dominated by
the elasticity of the bubbles, while at larger bubble sizes the foam
elasticity depends strongly on $\phi$.

As the foams coarsen the average bubble size increases and elastic effects
from the emulsion become relatively more important. In order to compare the
latter to Laplace pressure, we can construct a second elastocapillary number
as $Ca_{\text{el}}^{\text{FE}}=G_{\text{FE}}/(\gamma/R)=f(\epsilon)+\langle
R\rangle G_{0}\epsilon^{2}/\gamma$, which increases with $\langle R\rangle$
(Fig. 3d). It has been shown that monodisperse foams stop coarsening once
$Ca_{\text{el}}^{\text{FE}}\approx 0.8$ [9]. We can describe the evolution of
$R_{\text{sk}}$ with $G_{0}$ assuming a critical
$Ca_{\text{el}}^{\text{FE}}=0.54$ (Fig. 3c) unique to all $\phi$. Therefore,
in our foams elastic effects from the continuous phase become visible in the
bubble size distributions once $2G_{\text{FE}}\approx\gamma/\langle R\rangle$.
The experiment at $\phi=$ 0.65 stops at $Ca_{\text{el}}^{\text{FE}}\approx
0.4$, below the critical $Ca_{\text{el}}^{\text{FE}}$, which is why for these
foams the distribution are unchanged (Fig. 3a). The foams also never stop
coarsening despite reaching $Ca_{\text{el}}^{\text{FE}}=0.8$, however this is
because they are polydisperse.

Bubble size polydispersity in our samples also plays a crucial role in their
structural evolution. Since $G_{\text{FE}}$ decreases with $R$ (Fig. 3d), in a
polydisperse foam larger bubbles are easier to deform than smaller ones. This
spatial heterogeneity in the foam mechanical properties affects the coarsening
process causing a peculiar evolution of the foam structure. Bubbles will grow
preferentially towards weaker regions, thus towards other large bubbles. The
resulting heterogeneous ageing, which leads to a segregation between smaller
and larger bubbles, becomes visible after 96 hours of ageing (Fig. 4a) for
samples with $\phi\geq 0.75$. The clusters of small bubbles enclosed by the
chains of larger ones are striking at $\phi$ = 0.8 (Fig. 4b.)

The emulsion elasticity also changes the structure more locally. We observe
the emulsion bulging into the bubbles at the vertices (Fig. 4c). Such
inversion of curvature is another sign of the easier deformability of the
bubbles instead of the emulsion, and departure from capillary mediated
Plateau’s laws.

The inhomogeneous evolution of the foam due to the heterogeneity in mechanical
properties with bubble size resembles the elastic ripening recently observed
by Rosowski et al. at a microscopic scale [25]. In their work ripening is
arrested in elastic regions, and continues in the less elastic ones, where the
inhomogeneity is created by changing the polymer matrix while in our case the
inhomogeneity is inherently due to the variation in the response of the foam
with bubble size.

If the foams are left to evolve for a very long time, the larger bubbles will
keep growing until they start to coalesce, leaving the islands of small
bubbles behind (Fig. 4d). The photograph is reminiscent of structures observed
during the final stages of viscoelastic phase separation identified by
Tanaka[26] in systems with high dynamic asymmetry. He proposed its occurrence
in foams despite the difference in the coarsening processes from the polymer
phase separation or colloidal gelation previously studied, and the images
could be a first indication of such a process at millimetric length scales.

Figure 4: Foamed emulsion structure. a, Sample structure after 96 hours at
different $\phi$. From top left to bottom left, clockwise: $\phi$=0.70, 0.75,
0.80, 0.85. The edge of each frame is 4 cm. The corresponding skew of the
distributions is indicated by a star in Fig. 3b. At higher $\phi$ the foam
structure becomes highly unconventional with unrelaxed bubble shapes and
uneven emulsion distribution. b, Sample at $\phi$=0.80 after 5 days. Chains of
large bubbles and regions of small bubbles become clearly visible. Edge size
12 cm. c, Enlargement showing the emulsion bulge at the ends of the thin
Plateau borders between large bubbles. Edge size 3 cm. d, Sample having
$\phi$=0.80 made with sunflower oil, after 12 days. The edge of the frame is
13 cm. Segregation into regions of large and small bubbles is complete.

We show that viscoelasticity of the foam continuous phase leads to a slowing
down of coarsening and a radical change in foam structure. Once a critical
$Ca_{\text{el}}^{\text{FE}}$ is reached, the bubble size distributions show
the accumulation of smaller bubbles, which eventually segregate from the
larger ones. This is caused by the interplay between bubble size and foam
elasticity, which combined with bubble size polydispersity results in highly
heterogeneous materials. Stress build up in the unrelaxed structures will be
transferred to the foams if dried or solidified resulting in mechanical
weakness. The deviation from the traditional cellular foam structure can have
other implications on the material properties, so understanding and control of
the structural evolution is crucial for the safe design of foamy materials.

## MATERIALS AND METHODS

### Materials

Air bubbles and oil droplets (rapeseed oil from Brassica Rapa, Sigma-Aldrich)
are dispersed in an aqueous solution containing 30 g/L of sodium dodecyl
sulfate (SDS, from Sigma-Aldrich).

### Emulsion generation

Concentrated O/W emulsions are generated by mechanically mixing the oil and
the aqueous surfactant solution with the double syringe method. This technique
allows to control the oil volume fraction and ensures a good reproducibility
of the samples. The resulting droplet size distribution is measured with laser
diffraction with a Mastersizer 3000E (Malvern Panalytical) equipped with a
Hydro SM wet dispersion unit. The size distribution is polydisperse, with a
span, defined as $d(90\%)-d(10\%)/d(50\%)$, ranging from 1.2 to 1. The
surface-weighted mean drop diameter ranges between 3 and 5 $\mu$m. Both the
average size and the distribution width decrease with increasing $\phi$,
consistently with emulsification in turbulent flow [27].

### Emulsion rheology

The mechanical properties of concentrated emulsions are measured with a MCR302
rheometer (Anton Paar). The emulsion storage and loss moduli are obtained by
performing oscillatory strain sweep tests (frequency 1 Hz, amplitude from 10-5
to 1) in a cylindrical Couette geometry (CC27), with a gap of 1.1 mm. All
measurements are performed at (20.3 $\pm$ 0.1)∘C.

### Foam generation

Emulsions are foamed with a planetary kitchen mixer (Kenwood MultiOne 1000 W).
The mixing speed is gradually increased from the minimum to the maximum level,
and after 25 minutes of whipping, up to roughly 90% of air is incorporated
inside the emulsion, without affecting the droplet size distribution. The
final foam liquid fraction is measured by weight: a glass container of known
volume $V_{\text{foam}}$ is filled with freshly-made foamed emulsion and then
weighed, so that, if the gas density is assumed to be zero, the mass of the
foam liquid content $m_{\text{em}}$ is measured. We assume the emulsion
density to be given by the weighted sum of the density of its two components,
namely
$\rho_{\text{em}}\simeq\phi\rho_{\text{oil}}+(1-\phi)\rho_{\text{water}}$, so
that the foam liquid fraction is retrieved through the relation
$\varepsilon=m_{\text{em}}/(\rho_{\text{em}}V_{\text{foam}})$. The final foam
liquid fraction ranges between 9% and 13%, and slightly increases with $\phi$,
consistently with the increase of emulsion bulk viscosity [28].

### Imaging and image treatment

Freshly-generated foamed emulsions are carefully sandwiched between two square
glass plates (edge length 24 cm), separated by rubber joint of thickness 1 mm
which sets the cell gap. The cell is then placed between two square metal
frames which are screwed to keep the cell closed. The cell is then put under a
digital camera (Basler acA3800 - 14um, resolution 3840x2748 pixel) equipped
with a lens (TAMRON 16mm f/1.4), and a square of LED lights provides rather
uniform illumination from above. A photo of the foamed emulsion is taken every
3 minutes at early stage, and then every 30 minutes at late coarsening stage.

Custom-made MATLAB scripts are used for processing the images. A first image
pre-treatment is carried out by cropping the raw frames around a region of
interest, adjusting the contrast and obtaining the 2D foam skeleton through a
watershed algorithm. The skeletonised frames are then processed with a second
MATLAB script based on the built-in function regionprops to retrieve the area
of each gas cell, from which the equivalent radius is calculated as
$R=\sqrt{A/\pi}$.

## References

  * [1] J. Von Neumann.  Metal Interfaces (American Society for Metals, Cleveland, 1952),.
  * [2] Andrew M. Kraynik, Stephan A. Koehler, Howard A. Stone, and Sascha Hilgenfeldt.  An accurate von Neumann’s law for three-dimensional foams.  Physical Review Letters, 86(12):2685–2688, 2001.
  * [3] W. W. Mullins.  The statistical self-similarity hypothesis in grain growth and particle coarsening.  Journal of Applied Physics, 59(4):1341–1349, 1986.
  * [4] I.M. Lifshitz and V.V. Slyozov.  The kinetics of precipitation from supersaturated solid solutions.  Journal of Physics and Chemistry of Solids, 19(1):35–50, 1961.
  * [5] Zenaida Briceño-Ahumada, Amir Maldonado, Marianne Impéror-Clerc, and Dominique Langevin.  On the stability of foams made with surfactant bilayer phases.  Soft matter, pages 1459–1467, 2016.
  * [6] Mark Buchanan.  Liquid Crystal Foams: Formation and Coarsening.  arXiv, page 0206477, 2002.
  * [7] W Kloek, T van Vliet, and M Meinders.  Effect of bulk and interfacial rheological properties on bubble dissolution.  Journal of Colloid and Interface Science, 237(2):158–166, 2001\.
  * [8] A J Webster and M E Cates.  Osmotic Stabilization of Concentrated Emulsions and Foams.  Langmuir, 62(17):595–608, 2001.
  * [9] Houda Bey, Frédéric Wintzenrieth, Olivier Ronsin, Reinhard Höhler, and Sylvie Cohen-Addad.  Stabilization of foams by the combined effects of an insoluble gas species and gelation.  Soft Matter, 13(38):6816–6830, 2017.
  * [10] I. Lesov, S. Tcholakova, and N. Denkov.  Factors controlling the formation and stability of foams used as precursors of porous materials.  Journal of Colloid and Interface Science, 426:9–21, 2014.
  * [11] M Kogan, L Ducloue, J Goyon, X Chateau, O Pitois, and G Ovarlez.  Mixtures of foam and paste: suspensions of bubbles in yield stress fluids.  Rheologica Acta, 52(3):237–253, 2013.
  * [12] Lucie Ducloué, Olivier Pitois, Julie Goyon, Xavier Chateau, and Guillaume Ovarlez.  Coupling of elasticity to capillarity in soft aerated materials.  Soft Matter, 10(28):5093, 2014.
  * [13] Robert W. Style, Rostislav Boltyanskiy, Benjamin Allen, Katharine E. Jensen, Henry P. Foote, John S. Wettlaufer, and Eric R. Dufresne.  Stiffening solids with liquid inclusions.  Nature Physics, 11(1):82–87, 2015.
  * [14] F Gorlier, Y Khidas, and O Pitois.  Journal of Colloid and Interface Science Coupled elasticity in soft solid foams.  Journal of Colloid and Interface Science, 501:103–111, 2017.
  * [15] Alesya Mikhailovskaya, Véronique Trappe, and Anniina Salonen.  Colloidal gelation, a means to study elasto-capillarity effects in foam.  Soft Matter, pages 2249–2255, 2020.
  * [16] L.J. Gibson and M.F. Ashby.  Cellular solids: structure and properties.  Cambridge edition, 1997.
  * [17] T. G. Mason, J. Bibette, and D. A. Weitz.  Elasticity of compressed emulsions.  Physical Review Letters, 1995.
  * [18] C. D. Schimming and D. J. Durian.  Border-crossing model for the diffusive coarsening of two-dimensional and quasi-two-dimensional wet foams.  Phys. Rev. E, 96:032805, Sep 2017.
  * [19] Emilie Forel, Benjamin Dollet, Dominique Langevin, and Emmanuelle Rio.  Coalescence in two-dimensional foams: A purely statistical process dependent on film area.  Phys. Rev. Lett., 122:088002, Feb 2019.
  * [20] Alexandre Kabla and Georges Debrégeas.  Local stress relaxation and shear banding in a dry foam under shear.  Physical review letters, 90(25):258303, 2003.
  * [21] C Gay, P Rognon, D Reinelt, and F Molino.  Rapid Plateau border size variations expected in three simple experiments on 2D liquid foams.  European Physical Journal E, 34(1), 2011.
  * [22] James A. Glazier, Michael P. Anderson, and Gary S. Grest.  Coarsening in the two-dimensional soap froth and the large- Q potts model: A detailed comparison.  Philosophical Magazine B: Physics of Condensed Matter; Statistical Mechanics, Electronic, Optical and Magnetic Properties, 62(6):615–647, 1990.
  * [23] James A. Glazier, Steven P. Gross, and Joel Stavans.  Dynamics of two-dimensional soap froths.  Physical Review A, 36(1):306–312, 1987.
  * [24] Sylvie Cohen-Addad, Reinhard Höhler, and Olivier Pitois.  Flow in Foams and Flowing Foams.  Annual Review of Fluid Mechanics, 45(1):241–267, 2013.
  * [25] Kathryn A. Rosowski, Tianqi Sai, Estefania Vidal-Henriquez, David Zwicker, Robert W. Style, and Eric R. Dufresne.  Elastic ripening and inhibition of liquid–liquid phase separation.  Nature Physics, 16(4):422–425, 2020.
  * [26] Hajime Tanaka.  Viscoelastic phase separation.  Journal of Physics: Condensed Matter, 12:R207–R264, 2000.
  * [27] Slavka Tcholakova, Ivan Lesov, Konstantin Golemanov, Nikolai D. Denkov, Sonja Judat, Robert Engel, and Thomas Danner.  Efficient emulsification of viscous oils at high drop volume fraction.  Langmuir, 27(24):14783–14796, 2011.
  * [28] Nadya Politova, Slavka Tcholakova, Zhulieta Valkova, Konstantin Golemanov, and Nikolai D. Denkov.  Self-regulation of foam volume and bubble size during foaming via shear mixing.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 539(December 2017):18–28, 2018.

Acknowledgements
This study was partly supported by the Centre National d’Études Spatiale
(CNES) through project "Hydrodynamics of wet foams". A.S. acknowledges
Véronique Trappe and Eric Dufresne for interesting discussions.

Author contributions
C.G., E.R. and A.S. designed the experimental study. C.G. carried out the
experiments. C.G., J.M.I., E.R., A.P. and A.S. analysed and interpreted the
experimental results. C.G., E.R. and A.S. wrote the manuscript with
contributions from J.M.I and A.P.

Competing interests
The authors declare no competing financial interests.