CWCA: Complex-valued encoding water cycle algorithm

1.   Introduction

Generally, the objective of optimization is to find the optimal feasible response, considering the constraints of a problem. As a modern optimization method, meta-heuristic algorithms have been proposed by researchers in recent years, for example, the genetic algorithm (GA) ^[1], particle swarm optimization (PSO) ^[2], simulated annealing (SA) ^[3], glowworm swarm optimization (GSO) ^[4], harmony search (HS) ^[5], bacterial foraging optimization ^[6], earthworm optimization ^[7], bat algorithm (BA) ^[8], elephant herding behavior algorithm ^[9], krill herd algorithm ^[10], water cycle algorithm (WCA) ^[10], and various improved versions of the meta-heuristic optimization algorithm ^{[45,46,47,48]}. To date, hundreds of meta-heuristics have been proposed and used successfully to solve complex optimization problems.

In 2012, Eskandar et al. presented the WCA and used it to solve real-valued optimization issues. The real-valued WCA is a meta-heuristic optimization algorithm inspired by the water cycle process in nature. It considers how rivers and streams flow into the sea. A simplified water cycle system diagram is illustrated in Figure 1. In the WCA, the initial population comprises raindrops; the best raindrop represents the ocean. Many good raindrops represent a river, and the remaining raindrops represent streams that flow to the ocean and rivers. The rivers flow to the ocean, which is the lowest terrain location.

Since the WCA was proposed, it has been used extensively in scientific computation and engineering optimization. In 2014, Ail et al. applied the WCA to multi-objective optimization ^[12]. Zhang et al. applied the WCA to solve engineering optimization problems ^[13]. In 2015, Ail et al. proposed the WCA with the evaporation rate (ERWCA) for unconstrained and constrained optimization ^[14]. In ^[41], a comprehensive and exhaustive review was conducted on the WCA and its applications in a wide variety of study fields, including mechanical engineering, electrical and electronic engineering, civil engineering, industrial engineering, water resources and hydropower engineering, computer engineering, and mathematics. In ^[42], the detailed open source code for the WCA was provided, and its performance and efficiency for solving optimization problems was demonstrated. In ^[43], an enhanced discrete version of the WCA called DWCA was proposed to solve the symmetric and asymmetric traveling salesman problem. The designed solver was tested on over 33 problem datasets, and the statistical significance of the performance gaps for this benchmark was validated using results from non-parametric tests, not only in terms of optimality but also in terms of convergence speed. In ^[44], an extended version of WCA, that is, gradient-based WCA (GWCA) with the evaporation rate, was introduced to enhance the performance of the standard WCA by incorporating a local optimization operator in a so-called gradient-based approach. The experimental results demonstrated the feasibility and efficiency of the proposed GWCA.

In this paper, the complex-valued encoding WCA (CWCA) is proposed, which uses the complex number coding method of the complex-valued BA ^[15] and individual genes in the evolutionary algorithm ^[16,17]. A diploid is adopted in the indication of individual genes, which significantly enhances the individual's in information capacity. Finally, twelve benchmark functions and four engineering examples were considered to evaluate the performance of the CWCA. The test results indicated that the CWCA had higher precision and convergence speed than the real-valued WCA and other well-known meta-heuristic algorithms.

2.   Mathematical model for the real-valued water cycle

In the WCA, the initial population is represented by $N_{p o p} \times N$ matrix:

where ${N_{pop}}$ represents the population size and N represents the design variables. Initially, the ${N_{pop}}$ stream is randomly generated, and ${N_{sr}}$ good individuals represent a sea or river. The stream that has the minimum value among all streams is identified as the ocean.${N_{sr}}$ represents the summation of the number of rivers.

The remainder of the populations are computation using

The total volume of water flowing in a river and/or the ocean varies from stream to stream. The number of specified streams for a river or the ocean can be calculated as

where $N S_{n}$ denotes the number of rivers and f represents the fitness function.

In nature, streams are formed by raindrops, and then they join each other to constitute new rivers. A part of a streams flows directly into the ocean. All streams and rivers end up in the open ocean (the best point). Figure 2 shows a schematic of a stream flowing into a river. The star represents the river and the circles represent the stream.

In Figure 2, a stream flows into a river using a random distance:

where C is a numerical value between 1 and 2 (near to 2). The best numerical value is chosen as 2. d represents the current distance between the stream and river current distance. X is a random number between 0 and ($C \times d$).

The position update for streams and rivers is determined

where rand represents a random number uniformly distributed between 0 and 1. If the current solution provided by a stream is better than that of its connecting river, the positions of the river and stream are swapped. Such a swap can also occur for rivers and the ocean.

New raindrops form streams in different places to avoid being trapped in local optima. The rainfall condition is as follows:

where eps is a small numerical constant.

The stream and river location update formulas are

where LB represents the lower bound and UB represents the upper bound. The best newly formed raindrop is treated as a river running into the ocean. The remaining raindrops form new streams that flow into the rivers or ocean.

In Figure 3, the circles, stars, and diamonds represent streams, rivers, and the ocean, respectively. In view to the above description, the WCA steps are as follows:

Step 1. Initialize parameters ${N_{pop}}, {N_{sr}}$, and the number of maximum iterations iterMax

Step 2. Randomly generate the initial population, composite initial rivers, streams, and ocean using Eqs (1)–(3).

Step 3. Calculate the intensity of streams flowing into the rivers and ocean using Eq (4).

Step 4. Determine the streams' flow into the rivers, and the rivers' flow into the ocean using Eqs (6) and (7).

Step 5. According to the fitness function value, swap the positions of the streams, rivers, and ocean.

Step 6. If the termination condition is not satisfied [using Eq (8)], the process of raining occurs using Eqs (9) and (10); otherwise, go to Step 4.

Step 7. Stop the algorithm if the termination condition is satisfied; otherwise, go to Step 4.

3.   Complex-valued encoding WCA

3.1.   Complex-valued encoding method

In the CWCA, the imaginary and real parts of the complex number are updated respectively for each individual, which leads to inherent parallelism, improves the diversity of individuals, and enhances its exploitation and exploration capabilities, and it does not fall into the local optimum.

3.1.1.   Initialize the CWCA population

In the WCA, first, an interval $\left[L B_{k}, U B_{k}\right], k = 1, 2, \ldots, N_{p o p}$ is defined, and the ${N_{pop}}$ complex modulus and phase angle are randomly generated using

The 2M complex number is obtained as follows:

Therefore, in the CWCA, the initial population represented by a matrix of streams of size $N_{p o p} \times N$ is composed as follows:

where ${N_{pop}}$ denotes the complex modulus.

3.1.2.   Updating method of the CWCA

In the CWCA, the new positions of rivers and streams are determined as follows: the real parts are updated using

and the imaginary parts are updated using

3.1.3.   Fitness value

The fitness value is determined as follows: First, the complex value is converted into a real value and then its fitness value is calculated. The real-valued fitness function is determined using the complex modules and the sign is updated using the amplitude angle:

where $X_{k}$ denotes the converted real variables.

3.2.   Pseudocode of the CWCA

The pseudocode of the CWCA is as follows:

4.   Experimental results and comparisons

The experimental setup was as follows: MATLAB R2012a, AMD Athlon™*4 640 processor, and 2 GB memory. The CWCA was used to solve 12 benchmark functions ^[18,19] and four engineering design problems. The performance of CWCA was compared with that of GSO ^[4], artificial bee colony (ABC) ^[20], WCA ^[11], and ERWCA ^[14] using the standard deviation and mean. The control parameters of the algorithms were set as follows:

● GSO: Parameters ${N_{pop}}$ = 50, $G_{0} = 100, \alpha = 20; K_{0}$ linearly decreased to 1 and was set to NP ^[4].

● ABC: Parameters ${N_{pop}}$ = 50, Limit = 5D. ^[20].

● WCA: Parameters ${N_{pop}}$ = s 50, ${N_{sr}}$ = 8 (number of rivers), with C = 2. ^[11].

● ERWCA: Parameters ${N_{pop}}$ = 50, $N_{s r} = 8$, with C = 2. ^[14].

● CWCA: Parameters ${N_{pop}}$ = 50, $N_{s r} = 8$, with C = 2, $\rho_{k} \in\left[0, \frac{L B_{k}-U B_{k}}{2}\right], \theta_{k} \in[-2 \pi, 2 \pi]$

For the simulation experiments, the mean outcomes of 25 independent runs for f₀₁ - f₀₄ are listed in Table 2, the test results for functions f₀₅ - f₀₈ are listed in Table 3, and the test results for functionsf₀₉ - f₁₂ are listed in Table 4. The unrestrained benchmark functions are presented in Table 1. The results perform 25 generations. "Best, " "Mean, " "Worst, " "Std, " and "ANOVA" represent the optimal value, mean value, worst value, standard deviation, and analysis of variance, respectively.

4.1.   Experimental results for the benchmark functions

As shown in Table 2, the best values produced by the CWCA were more accurate than those obtained using the proposed WCA and the other meta-heuristic optimization algorithms. For four functions, the standard deviation of the CWCA was less than that of GSO, ABC, WCA, and ERWCA, which implies that the CWCA had better stability in terms of optimizing high-dimensional unimodal functions. Figures 4–7 show that the fitness functions converged to a curve. Hence, it can be concluded that the proposed CWCA had a higher convergence rate and higher computation precision than meta-heuristic algorithms such as GSO, ABC, WCA, ERWCA, and CWCA. Figures 8–11 show the analysis of ANOVA using a graph. The ERWCA and CWCA obtained stable optimal values. Figure 9 show that most of the algorithms produced a stable value, except the WCA, when solving function f₀₂.

It can be noted from Table 3 that the CWCA located the optimal solutions of the functions f₀₅, f₀₆ and f₀₈ with a standard deviation of 0. For function f₀₇, the precision and mean value were higher than those for the other meta-heuristic optimization algorithms. Building on the results shown in Figures 12-15, it can be concluded that the CWCA converged faster and its precision was higher than that of the other meta-heuristic optimization algorithms. Figures 16-19 show the ANOVA graphical analysis results. The results demonstrated that both the ERWCA and CWCA had better stability than the other meta-heuristic optimization algorithms.

As shown in Table 4, for f₀₉, all algorithms produced optimal solutions, and the solutions of GSO were the most stable. For functions f₁₀ and f₁₁, the solutions for most algorithms achieved the same order of magnitude, but the standard deviation of the CWCA was the best. For function f₁₂, both the ERWCA and CWCA obtained optimal solutions, but the CWCA had the best stability. Figures 20–27 shows that the CWCA had excellent performance when solving multimodal low-dimensional problems.

4.2.   Engineering design problem

The mathematical model for general constrained optimization is defined ^[21] as

where $x = \left(x_{1}, x_{2}, \ldots, x_{D}\right)$ represents the D-dimensional decision variables, $f(x){\rm{ }}$ denotes the objective function, ${g_j}(x) \le 0$ represents q inequality constraints, and ${h_i}(x)$ represents m - q equality constraints. The functions f, g_j and ${h_i}$ are nonlinear or linear functions. ${l_i}, {u_i}$ denote the lower and upper bounds, respectively.

Generally, when practical problems are solved, it is common to convert equality constraints into inequality constraints:

where $\varepsilon$ denotes the allowed tolerance, which is set to a small value. The best numerical value for $\varepsilon$ is 0.001. According to the feasibility-based rules, each individual's constraint violation is

Generally, a feasibility-based rule is used ^[22]:

Rule 1: Any feasible solution is preferred over any infeasible solution.

Rule 2: The feasible solution that has a better objective function value is preferred, if choosing between two feasible solutions.

Rule 3: The infeasible solution that has smaller constraint violations is preferred, if choosing between two infeasible solutions.

4.2.1.   Pressure vessel design problem

A cylindrical pressure vessel with hemispherical heads capped at both ends is shown in Figure 28, where $T_{s}\left(x_{1}\right)$ denotes the shell thickness, $T_{h}\left(x_{2}\right)$ denotes the head thickness, $R\left({{x_3}} \right)$ denotes the inner radius, and $L\left({{x_4}} \right)$ denotes the cylindrical vessel length. The overall cost, including a combination of material cost, single 60⁰ welding cost, and forming cost, is minimized.

The mathematical model for the pressure vessel design problem is

The results obtained using the CWCA to solve the pressure vessel design problem is shown in Table 5. The statistical outcomes were compared with those of the meta-heuristic algorithms ES ^[27], GSA ^[28], PSO ^[29], GA ^[24,25,26], DE ^[30], MVO ^[23], WCA ^[11], and ACO ^[31]. As shown in Table 5, the proposed CWCA outperformed all the other optimization algorithms.

4.2.2.   Cantilever beam design problem

Cantilever beam design optimization aims to minimize the weight of a cantilever beam, which is composed of hollow square blocks. In this study, it involved five squares: the first block was fixed, and the fifth block had the burden of a vertical load. In Figure 29, five parameters define the shape of the cross section of the cubes.

The mathematical model for the cantilever beam design optimization is

In the same way, the meta-heuristic optimization algorithms MVO ^[23], CS ^[32], SOS ^[33], MMA ^[34], WCA ^[11], and CWCA were used. Table 6 shows that the results of the CWCA for the cantilever beam design problem were consistent with those of the other engineering optimization problems. Both the WCA and CWCA outperformed all the other optimization algorithms.

4.2.3.   Three-bar truss design problem

The design problem for a three-bar planar truss structure is shown in Figure 30. The volume of a statically loaded three-bar truss was minimized, subject to stress (σ) constraints. The objective was to find the best cross-sectional areas (A₁, A₂).

The mathematical model for the three-bar planar truss structure design problem is

For the three-bar truss design problem, the metaheuristic optimization algorithms MVO ^[23], DEDS ^[34], PSO-DE ^[35], MBA ^[36], CS ^[37], WCA ^[11], and CWCA were used. The proposed CWCA was compared with the other optimization algorithms (see Table 7). The results obtained using the CWCA were very close to those obtained using DEDS.

4.2.4.   Welded beam design problem

Figure 31 shows a beam made of low-carbon steel and welded to a rigid support. The welded beam was designed to achieve the minimum cost, subject to constraints on shear stress (τ), bending stress (σ) in the beam, buckling load on the bar (P_b), end deflection of the beam (δ), and side constraints. There were four design variables: $h\left({{x_1}} \right), l\left({{x_2}} \right), t\left({{x_3}} \right)$ and $b\left({{x_4}} \right)$.

The objective function is expressed as

where

In the same way, the meta-heuristic optimization algorithms MVO ^[23], GSA ^[28], CPSO ^[38], GA ^[39], and HS ^[40] were applied. The proposed CWCA was compared with the other meta-heuristic algorithms (see Table 8). The results showed that the proposed CWCA achieved the minimum cost.

5.   Conclusions and future directions

In recent years, significant attention has been paid to the design of meta-heuristic optimization algorithms to solve optimization problems. The aims of this study were to improve the exploration and exploitation capabilities of the WCA algorithm, accelerate its convergence speed, and enhance its calculation accuracy. A novel CWCA algorithm was proposed, and to evaluate the performance of CWCA optimization, 12 benchmark functions and four engineering examples were considered. The test results indicated that the CWCA had higher precision and convergence speed than the WCA and other popular meta-heuristic optimization algorithms. Further improvements to versions of the WCA equipped with the latest, efficient strategies should be considered in future research. Additionally, the CWCA is not well used in some engineering majors, such as mechanical engineering. Despite the great potential of this principal field, it is necessary to use and use more meta-heuristic optimization algorithms, such as CWCA, to solve large-scale optimization problems.

Acknowledgments

This work was supported by the National Science Foundation of China under Grant No. 62066005 and the Project of Guangxi Natural Science Foundation under Grant No. 2018GXNSFAA138146. We thank Maxine Garcia, Ph.D, from Liwen Bianji (Edanz) (www.liwenbianji.cn/) for editing the English text of a draft of this manuscript.

Conflict of interest

The authors declare no conflict of interest.

Appendix

CWCA MATLAB source code

%-----------------------------------------------

function [Xmin, Fmin] = CWCA(objective_function, LB, UB, nvars, Fnum, Nmax)

  nvar = nvars;    % Dimension

  N = 50;    % Population size

  Nsr = 8;    % Number of rivers

  C = 2;

  Nwat = N-Nsr;    % Raindrops

  Fmin = inf;

%------------- Initialization (In module and Angle mode)----------------------------------

  zigma_NCn = inf*ones(1, Nwat + 1);

for i = 1: N

  x_age(i, :) = 4*pi*(rand(1, nvar)-0.5); % Argument

  x_mod(i, :) = 0.5*(UB - LB).*rand(1, nvar); % Module

  x_R(i, :) = x_mod(i, :).*cos(x_age(i, :));

  x_I(i, :) = x_mod(i, :).*sin(x_age(i, :));

  sq_x(i, :) = sqrt(x_R(i, :).^2 + x_I(i, :).^2);

  x(i, :) = sign(sin(x_I(i, :)/sq_x(i, :)))*sq_x(i, :) + 0.5*(UB + LB);

end

  LB_bound = repmat(LB, N, 1);

  UB_bound = repmat(UB, N, 1);

  x = max(x, LB_bound);

  x = min(x, UB_bound);

for i = 1 : N

  xx(i, :) = objective_function(x(i, :), Fnum);

end

  [~, index1] = sort(xx, 'ascend');

  Sr = x(index1(1: Nsr+1), :);

  sr_R = x_R(index1(1: Nsr+1), :);

  sr_I = x_I(index1(1: Nsr+1), :);

  SR = sr;

  SR_R = sr_R;

  SR_I = sr_I;

%--------------------------------------------------------------------------

for i = 1:Nsr+1

  cost(i) = objective_function(SR(i, :), Fnum); % Calculate fitness values

end

  cs = sort(cost, 'ascend')';

  cs = cs(1: Nsr+1);

  CN = cs-max(cs);

  Pn = abs(CN/(sum(CN)));

  Pn(Nsr+1) = [];

  NCn = round(Nwat*Pn);

while (sum(NCn) ~ = Nwat)

  i = Nsr;

while sum(NCn) > Nwat

  if NCn(i) < = 1

    i = i-1;

end

    NCn(i) = NCn(i)-1;

end

    i = 1;

while sum(NCn) < Nwat

    NCn(i) = NCn(i)+1;

end

end

  NCn = sort(NCn, 'descend');

  Sr = sr(1: Nsr, :);

  sr_R = sr_R(1: Nsr, :);

  sr_I = sr_I(1: Nsr, :);

  zigma_NCn = [0];

for i = 1:Nsr

  zigma_NCn = [zigma_NCn sum(NCn(1: i))];

end

end

%------------- Initialization ----------------------------------

for j = 1:Nwat

  water_age(j, :) = 4*pi*(rand(1, nvar)-0.5);

  water_mod(j, :) = 0.5*(UB - LB).*rand(1, nvar);

  water_R(j, :) = water_mod(j, :).*cos(water_age(j, :));

  water_I(j, :) = water_mod(j, :).*sin(water_age(j, :));

  sq_water = sqrt(water_R(j, :).^2 + water_I(j, :).^2);

  water(j, :) = sign(sin(water_I(j, :)/sq_water))*sq_water + 0.5*(UB + LB);

end

  WATER_R = water_R;

  WATER_I = water_I;

  WATER = water;

%------------------------------------------------------

  Sea = sr(1, :);

  sea_R = sr_R(1, :);

  sea_I = sr_I(1, :);

%--------------- The main loop ------------------------------

  Locate = 1;

for ii = 1:Nmax % Maximum number of iterations

%--------------- The stream flows to the river --------------------------

  stp = 1;

for i = 1: Nsr_k

  for j = (zigma_NCn(i) + 1): zigma_NCn(i+1)

%----------- Real component update --------

  new_WATER_R = WATER_R(j, :) + C.* rand(1, nvar).*(sr_R(i, :) - WATER_R(j, :));

%----------- Imaginary part update --------

  new_WATER_I = WATER_I(j, :) + C.* rand(1, nvar).*(sr_I(i, :) - WATER_I(j, :));

  sq_new_WATER = sqrt(new_WATER_R.^2 + new_WATER_I.^2);

  new_WATER =

sign(sin(new_WATER_I/sq_new_WATER))*sq_new_WATER+0.5*(UB+LB);

if objective_function(WATER(j, :), Fnum) > objective_function(new_WATER, Fnum)

  WATER(j, :) = new_WATER;

  WATER_R(j, :) = new_WATER_R;

  WATER_I(j, :) = new_WATER_I;

end

end

end

%------------Rivers flow to the sea --------------------------

for i = 1:Nsr_k

  new_sr = [];

%----------- Real component update --------

  new_sr_R = sr_R(i, :) + C.*rand(1, nvar).*(sea_R - sr_R(i, :));

%----------- Imaginary part update --------

  new_sr_I = sr_I(i, :) + C.*rand(1, nvar).*(sea_I - sr_I(i, :));

  sq_new_sr = sqrt(new_sr_R.^2 + new_sr_I.^2);

  new_sr = sign(sin(new_sr_I/sq_new_sr))*sq_new_sr + 0.5*(UB + LB);

if objective_function(sr(i, :), Fnum) > objective_function(new_sr, Fnum)

  sr_R(i, :) = new_sr_R;

  sr_I(i, :) = new_sr_I;

  sr(i, :) = new_sr;

 end

end

%------------ Replacement of rivers and streams --------------------

  new_sr = [];

for i = 1:Nsr_k

   for j = (zigma_NCn(i)+1):zigma_NCn(i+1)

    if objective_function(sr(i, :), Fnum) > objective_function(WATER(j, :), Fnum)

temp = WATER(j, :);

      temp_R = WATER_R(j, :);

      temp_I = WATER_I(j, :);

      WATER(j, :) = sr(i, :);

      WATER_R(j, :) = sr_R(i, :);

      WATER_I(j, :) = sr_I(i, :);

      sr(i, :) = temp;

      sr_R(i, :) = temp_R;

      sr_I(i, :) = temp_I;

      end

    end

  end

%------------- River and sea location update ----------------------

for i = 1:Nsr_k

   if objective_function(sea, Fnum) > objective_function(sr(i, :), Fnum)

      temp = sr(i, :);

      temp_R = sr_R(i, :);

      temp_I = sr_I(i, :);

      sr(i, :) = sea;

      sr_R(i, :) = sea_R;

      sr_I(i, :) = sea_I;

      sea = temp;

      sea_R = temp_R;

      sea_I = temp_I;

   end

end

   Xmin = sea;

   Fmin = objective_function(sea, Fnum);

   f(ii) = Fmin;

end

end

%----------------------End------------------------------------------------