Artificial Protozoa Optimizer (APO): A novel bio-inspired metaheuristic algorithm for engineering optimization

In plain words

Protozoa are tiny single-celled animals. This method copies a few of their habits to search for the best answer to a problem. Picture a crowd of candidate answers as a group of protozoa. Each one tries to move to a better spot. Some move toward a nearby spot, and others move by mixing in the positions of their neighbours in the ranking. The weaker candidates sometimes go dormant and are replaced by a random new position. The others sometimes split into a copy that is nudged a little. In the algorithm as printed, a candidate only moves if its new spot scores better.

The method has only two settings of its own. The authors tested it on a standard set of 12 benchmark functions, on five classic engineering design problems, and on splitting one colour image into regions by pixel value.

TL;DR

• APO is a population-based metaheuristic for black-box minimisation, with search moves modelled on protozoan foraging, dormancy and binary fission §2.2.
• It has two method-specific parameters: the number of neighbour pairs np and the maximum dormancy/reproduction fraction pf_max §2.2.4.
• Author-reported: ranked first by Friedman test among 17 algorithms on CEC2022 at dimension 20 §3.3.1 tab. 5, and first by mean rank over the five engineering designs tab. 11.
• Author-reported: first by mean rank on a single-image multilevel thresholding test (Lena, 2 to 10 thresholds) tab. 14.
• Use it as a baseline for black-box problems. The evidence is one benchmark suite, five engineering problems and one image; the code was not verified (see Paper vs code).

Problem

Exact optimisation methods need derivatives or convexity and can get stuck in local optima. Metaheuristics search without those assumptions but give no optimality guarantee §1, pp. 1–2. The authors frame APO as a new population method for this setting §1, p. 3.

Method

APO keeps ps candidate positions ("protozoa") in dim dimensions. Each iteration sorts the population by fitness, then updates every protozoan once. Each one either goes dormant or reproduces (a random subset chosen by pf), or forages. Foraging is either autotrophic or heterotrophic, chosen by a probability that decays over the run §2.2.4 alg. 1. A new position replaces the old one only if its fitness is lower (greedy acceptance, as printed in the algorithm).

Steps for one iteration iter (1-based), following alg. 1:

1. Sort the protozoa by fitness, smallest first alg. 1 l.2 tab. 1.
2. Draw pf = pf_max · rand eq. 15 and pick ceil(ps · pf) ranks at random for the dormancy/reproduction set alg. 1 l.3–4.
3. For a protozoan in that set: dormancy with probability p_dr eq. 17, giving a uniform random position in the bounds eq. 11; otherwise reproduction, a perturbation on a random subset of dimensions eq. 13 eq. 14.
4. For any other protozoan: build a foraging mask of ceil(dim · i / ps) dimensions by rank eq. 6 alg. 1 l.15–16. With probability p_ah eq. 16 use the autotroph move eq. 1, otherwise the heterotroph move eq. 7.
5. Accept the candidate if its fitness is lower than the current one alg. 1 l.23–27.
6. Record the best protozoan alg. 1 l.29 and count ps evaluations alg. 1 l.30.

The foraging factor f starts near 2 and decays to 0 over the run eq. 3 fig. 6.

# Artificial Protozoa Optimizer (APO)  [alg. 1]
# Input:  objective f; bounds X_min, X_max (dim,); ps, np, pf_max, dim,
#         iter_max
# Output: best position X_gbest and its fitness
# Ranks i = 1..ps, as in the paper. Minimisation. The helpers read ps,
# np, dim, X_min, X_max, iteration and iter_max from the enclosing scope.

def apo(f, X_min, X_max, ps, np, pf_max, dim, iter_max):
    # Random initial population, one evaluation each  [§2.3.1]
    X = X_min + rand(ps, dim) * (X_max - X_min)     # (ps, dim)
    fit = f(X)                                       # (ps,)
    for iteration in range(1, iter_max + 1):
        # Rank by fitness, smallest first  [alg. 1 l.2] [tab. 1]
        X, fit = sort_by_fitness(X, fit)
        # Share of dormant/reproducing protozoa  [eq. 15]
        pf = pf_max * rand()
        # Random ranks for dormancy/reproduction  [alg. 1 l.4]
        dr_index = randperm(ps, ceil(ps * pf))
        for i in range(1, ps + 1):
            # Foraging factor, decays from 2 to 0  [eq. 3]
            forage = rand() * (1 + cos(iteration / iter_max * pi))
            # Autotroph probability, decays from 1 to 0  [eq. 16]
            p_ah = 0.5 * (1 + cos(iteration / iter_max * pi))
            if i in dr_index:
                candidate = dormant_or_reproduce(X, i)
            else:
                candidate = forage_move(X, fit, i, forage, p_ah)
            # Greedy acceptance, same test for every branch  [alg. 1 l.23]
            value = f(candidate)
            if value < fit[i]:
                X[i] = candidate
                fit[i] = value
        # Best protozoan after this iteration  [alg. 1 l.29]
        X_gbest = X[argmin(fit)]
        # ps evaluations per iteration  [alg. 1 l.30]
    return X_gbest, min(fit)

def dormant_or_reproduce(X, i):
    # Dormancy probability grows with the rank i  [eq. 17]
    p_dr = 0.5 * (1 + cos((1 - i / ps) * pi))
    if p_dr > rand():
        # Dormant: a fresh random position in the bounds  [eq. 11]
        return X_min + rand(dim) * (X_max - X_min)
    # Reproduction: perturb X_i on ceil(dim * rand) dimensions  [eq. 14]
    mask_r = zeros(dim)
    mask_r[randperm(dim, ceil(dim * rand()))] = 1
    # Forward or reverse perturbation  [eq. 13]; inner term: see Pitfalls
    sign = random_choice([+1, -1])
    step = X_min + rand(dim) * (X_max - X_min)
    return X[i] + sign * rand() * step * mask_r

def forage_move(X, fit, i, forage, p_ah):
    # Foraging mask: ceil(dim * i / ps) dimensions  [eq. 6] [alg. 1 l.16]
    mask_f = zeros(dim)
    mask_f[randperm(dim, ceil(dim * i / ps))] = 1
    if p_ah > rand():
        # Autotroph: pull toward a random partner  [eq. 1] [§2.2.1]
        j = randint(1, ps)
        pull = zeros(dim)
        for k in range(1, np + 1):
            # Pair k: a lower rank and a higher rank, rank 1 and ps
            # use the protozoan itself  [§2.2.1]
            lower = randint(1, i - 1) if i > 1 else 1
            upper = randint(i + 1, ps) if i < ps else ps
            # Pair weight from the fitness ratio  [eq. 5]
            w_a = exp(-abs(fit[lower] / (fit[upper] + eps)))
            pull = pull + w_a * (X[lower] - X[upper])
        return X[i] + forage * (X[j] - X[i] + pull / np) * mask_f
    # Heterotroph: move around a nearby spot  [eq. 7] [eq. 8]
    sign = random_choice([+1, -1])                  # see Pitfalls
    X_near = (1 + sign * rand(dim) * (1 - iteration / iter_max)) * X[i]
    pull = zeros(dim)
    for k in range(1, np + 1):
        # Fixed ranks i-k and i+k, clipped to 1..ps  [§2.2.1]
        lower = clip(i - k, 1, ps)
        upper = clip(i + k, 1, ps)
        # Pair weight  [eq. 9]
        w_h = exp(-abs(fit[lower] / (fit[upper] + eps)))
        pull = pull + w_h * (X[lower] - X[upper])
    return X[i] + forage * (X_near - X[i] + pull / np) * mask_f

Key equations, as printed (the operator ⊙ is element-wise multiplication):

X_i^{new} = X_i + f \cdot \Big(X_j - X_i + \frac{1}{np} \sum_{k=1}^{np} w_a \cdot (X_{k-} - X_{k+})\Big) \odot M_f eq. 1

f = rand \cdot \Big(1 + \cos\big(\tfrac{iter}{iter_{max}} \cdot \pi\big)\Big) eq. 3

w_a = \exp\Big(-\Big|\frac{f(X_{k-})}{f(X_{k+}) + eps}\Big|\Big) eq. 5

M_f[d_i] = \begin{cases} 1 & \text{if } d_i \in randperm(dim, \lceil dim \cdot i / ps \rceil) \\ 0 & \text{otherwise} \end{cases} eq. 6

X_i^{new} = X_i + f \cdot \Big(X_{near} - X_i + \frac{1}{np} \sum_{k=1}^{np} w_h \cdot (X_{i-k} - X_{i+k})\Big) \odot M_f eq. 7

X_{near} = \Big(1 \pm Rand \cdot \big(1 - \tfrac{iter}{iter_{max}}\big)\Big) \odot X_i eq. 8

w_h = \exp\Big(-\Big|\frac{f(X_{i-k})}{f(X_{i+k}) + eps}\Big|\Big) eq. 9

X_i^{new} = X_{min} + Rand \odot (X_{max} - X_{min}) eq. 11

X_i^{new} = X_i \pm rand \cdot \Big(X_{min} + Rand \odot (X_{max} - X_{min})\Big) \odot M_r eq. 13 (inner term read from a low-resolution page image; see Pitfalls)

pf = pf_{max} \cdot rand eq. 15

p_{ah} = \tfrac{1}{2}\Big(1 + \cos\big(\tfrac{iter}{iter_{max}} \cdot \pi\big)\Big) eq. 16

p_{dr} = \tfrac{1}{2}\Big(1 + \cos\big((1 - \tfrac{i}{ps}) \cdot \pi\big)\Big) eq. 17

Symbols: rand is one uniform draw in [0,1]; Rand is a vector of dim such draws eq. 10; eps = 2.2204e−16 §2.2.1; np_max = ⌊(ps − 1)/2⌋ is only the largest allowed np eq. 4. Mapping vectors M_f and M_r are 0/1 masks of size (1 × dim). M_r selects ceil(dim · rand) random dimensions eq. 14.

Neighbour sampling (foraging). The autotroph pair X_{k-}, X_{k+} is random, drawn from ranks below and above i; if i is 1 or ps, the protozoan itself is used §2.2.1. The heterotroph pair uses fixed ranks i−k and i+k, clipped to the ends §2.2.1.

Bounds. X_min and X_max are the per-dimension bounds lb and ub eq. 12. Clipping after a move is not described.

Complexity. Sorting costs O(ps log ps). Foraging costs O((ps − ⌈ps·pf⌉) · np · dim), dormancy and reproduction cost O(⌈ps·pf⌉ · dim), and fitness evaluation costs O(ps · f(·)) per iteration. With np = 1 the total is O(iter_max · ps · (log ps + dim + f(·))) §2.3.4.

Exploration and exploitation diagnostic. Not part of the optimiser. The authors measure the population's spread div (mean absolute distance from the per-dimension median) and define the exploration rate err = div / div_max and the exploitation rate eir = (div_max − div) / div_max eq. 18 eq. 19 eq. 20 §3.3.1:

div = \frac{1}{dim \cdot ps} \sum_{j=1}^{dim} \sum_{i=1}^{ps} \big| median(x^j) - x_i^j \big| eq. 20

Hyperparameters

Results

Author-reported. Setup: CEC2022, 12 functions at dimension 20, 1,000,000 evaluations, 30 runs, compared with 16 algorithms in the main text (Table 3) and 16 more in the supplement §3.3 tab. 3. The benchmark functions:

Engineering success rates and run costs (Table 13). The authors report 100% success for APO on the spring, welded beam, speed reducer and three-bar problems, with 9,900, 11,255, 6,148 and 6,423 mean evaluations respectively tab. 13 §4.1.6. See Pitfalls for the pressure-vessel entry.

Feasible values reported in Table 12: spring 0.01273171; pressure vessel 6117.896; welded beam 1.765002; speed reducer 2996.868; three-bar truss 263.8963 tab. 12.

Figures (described, not reproduced): on most CEC functions the exploration rate falls and the exploitation rate rises over the run, with F7 and F8 as exceptions fig. 7 §3.3.1. In the two-dimensional plots the sampled points are sparse far from the optimum and dense near it fig. 8. The mean-fitness convergence curves show APO converging faster than the comparison methods on most functions fig. 9 §3.3.2.

Multilevel thresholding: the authors use the minimum cross-entropy objective §4.2.1 eq. 24 eq. 25 eq. 26 eq. 27 and one 512 × 512 colour image, with each channel segmented separately §4.2.2. Thresholds are rounded at each iteration. With 2 thresholds all nine compared algorithms tie on PSNR at 23.5625 tab. 14. The advantage appears at larger threshold counts.

C(P, Q) = \sum_{i=1}^{n} p_i \log \frac{p_i}{q_i} eq. 24

\{t_1^*, \dots, t_n^*\} = \arg\min f(t_1, \dots, t_n) eq. 25

f(t_1, \dots, t_n) = \sum_{i=1}^{L} i \cdot z(i) \cdot \log(i) - \sum_{k=1}^{n+1} \sum_{i=t_{k-1}}^{t_k - 1} i \cdot z(i) \cdot \log\big(u(t_{k-1}, t_k)\big) eq. 26

u(t_{k-1}, t_k) = \frac{\sum_{i=t_{k-1}}^{t_k - 1} i \cdot z(i)}{\sum_{i=t_{k-1}}^{t_k - 1} z(i)} eq. 27

Here t_0 = 1, t_{n+1} = L + 1, 0 < t_1 < … < t_n < L, z(i) is the number of pixels with value i, and L is the largest pixel value §4.2.1. The image metrics are PSNR, SSIM and FSIM §4.2.2.

Paper vs code

The official code could not be read, so this section records what was found and what was not checked. Details:

• The abstract points to the author's project page (seyedalimirjalili.com/projects) and to a MathWorks File Exchange entry. The project page lists a MATLAB zip of the code and states no licence. The File Exchange URL returned HTTP 404 on our fetch.
• The GitHub repository xiaopengzie/Metaheuristics (commit 9ff25f631e498dda772ad30fa68762319e66cb8f) has a README that cites this paper with its DOI, and says it holds the paper's MATLAB and Python code. The owner's GitHub profile has no name, so the link to the authors rests on the README alone. It holds Artificial Protozoa Optimizer/Artificial Protozoa Optimizer.zip. The repository has no top-level LICENSE file. The MAAPO/ folder has its own LICENSE, which belongs to a different method.
• We could not extract that zip (the shell extraction was denied), so no code file was read. No [code: …] citation is given, and the algorithm was not compared line by line against code.
• A GitHub repository search for "artificial protozoa optimizer" returned five repositories: Pouya-Azizzadh (Python port), ish59, Dingyihao-nuist (3D coverage application), nikameskhi (APO plus variants and a spring benchmark), and xiaopengzie. The first four were not used. None lists a licence.

Not checked, so the code-vs-paper differences are unknown: foraging and reproduction details, the dormancy replacement rule, the eq. 13 inner term, neighbour boundary handling, the iteration schedule, the constraint penalty, and the segmentation rounding.

Pitfalls

• Sorting key (eq. 2). The paper writes X_i = sort(X_i) without a key. Table 1 says ranking is by fitness, smallest first, and Algorithm 1 line 2 sorts the population that way. The ranks decide the mask size in eq. 6, the neighbour choice in eqs 1 and 7, and the p_dr value in eq. 17. Sort by fitness, not by coordinates tab. 1 alg. 1 l.2.
• Dormancy acceptance. The text says a dormant protozoan is replaced by its new position §2.2.2 fig. 4(c). Algorithm 1 applies the same greedy test to all branches, so a dormant candidate that is worse than the current position is rejected alg. 1 l.23–27. Follow the algorithm; the text reads as unconditional.
• Random versus ranked dormancy. Section 2.3.3 says inferior protozoa tend to go dormant and superior ones reproduce §2.3.3. Algorithm 1 draws the dormancy/reproduction set Dr_index uniformly at random, and p_dr uses the rank only inside that set alg. 1 l.4 eq. 17. The effective behaviour is closer to the algorithm than to the text.
• In-place updates. Algorithm 1 writes each X_i as it goes, so later protozoa read neighbours already updated in the same iteration. The paper does not say whether neighbours come from the start-of-iteration population alg. 1 l.5–28. The algorithm as printed implies in-place updates.
• Neighbour boundaries. Ranks 1 and ps reuse the protozoan itself for the missing neighbour §2.2.1. Use this exactly, or the first and last protozoa will behave differently.
• Division by fitness. Weights are exp(−|f_a / (f_b + eps)|) eq. 5 eq. 9. The benchmark optimum is shifted to 0 §3.3. When the fitness in the denominator is near 0, the ratio can be very large and the weight near 0. The paper does not discuss negative fitness or near-zero denominators beyond eps.
• Eq. 13 inner term. The reproduction equation is printed with a ± sign (the text confirms forward and reverse perturbations §2.2.3). The inner expression in this card was read from a low-resolution page image and is not confirmed. The PDF is the only source we have; check the inner term and the sign rule against the authors' files before implementing.
• Sign choice. Eq. 8 and eq. 13 use ±. The paper does not say whether the sign is drawn once per protozoan or per dimension §2.2.2 §2.2.3.
• Ceiling brackets. Eqs 6 and 14 use square brackets that Table 1 defines as ceiling. Use ceil: with a floor, dim · i / ps gives zero dimensions for small ranks, and the protozoan would not move during foraging tab. 1.
• Global best. Algorithm 1 line 29 uses opt(·), which is not defined. Take the minimum-fitness protozoan.
• Schedules depend on iter_max. Foraging factor and autotroph probability use iter/iter_max. The CEC runs use 10,000 iterations at ps = 100, but the engineering and segmentation runs use 500 and 100 iterations. The same code gives different schedules under each setting §3.3 §4.1 §4.2.2.
• Penalty coefficients. The engineering problems use a penalty function, but its coefficients are not stated §4.1. Results depend on them, so reproduce with the authors' code or state your own coefficients.
• Pressure-vessel entry in Table 13. The APO success rate for the pressure vessel is printed as 0.7419, without a % sign. The other four APO entries are 100%, and the text does not discuss this value. It may be a typo for 74.19, or it may be correct. Check the source before relying on it. Table 12's pressure-vessel feasible value (6117.896) is also not the Table 7 optimum (5887.614).
• Count of compared algorithms. The abstract says 32 algorithms. Table 5 and Table 3 each list 16 comparators, and the supplement holds 16 more §3.3. The supplement is not in our PDF, so only the 16 in the tables can be checked here.
• Run counts differ by section. CEC2022 uses 30 runs §3.1; engineering and segmentation use 31 §4.1 §4.2.2.
• One image for segmentation. The segmentation test is a single 512 × 512 colour image §4.2.2. Generalisation beyond that is not shown.
• Data availability. The authors state "data will be made available on request" p. 18. The code location is in Paper vs code.
• Cost. The complexity expression in §2.3.4 counts sorting, updates and fitness calls, not neighbour selection or the extra fitness lookups for weights. Wall-clock time versus other methods is not reported.

Lineage

• Builds on: no earlier method is named as its base. The comparison set is GA, DE, BSA, Jaya, PSO, GWO, WOA, PPE, GSA, MVO, SCA, AOA, TLBO, SDO, SPO and SPBO tab. 3.
• Compared against: the 16 main-text algorithms in Table 3 and the 16 in the supplement (Table A.3–A.6, not in our PDF) §3.3.
• Superseded by: not known from content in this wiki.