Pray Algorithm Optimization (PAO): An Astronomically Inspired Metaheuristic Based on Solar Altitude and Islamic Prayer Times for Global Optimization and Engineering Design

1. Introduction

This section explains why the paper exists. It describes what optimization problems are and why so many methods have been invented to solve them, then reviews the main families of these methods and how researchers group them. The review shows a gap: although algorithms have been drawn from animals, physics, chemistry, human society and mathematics, almost none has been built from the daily movement of the Sun. The section ends by stating what this paper contributes — a new optimizer called Pray Algorithm Optimization (PAO), built on the timing of the five daily prayers — and by listing the tests used to check whether it works.

1.1. Optimization and Its Reach

Optimization is, at its core, the search for the best available option among many. It appears wherever a resource is limited — time, money, raw material, energy, computing power — and wherever some quantity must be pushed as high or as low as the situation allows. Because limited resources are a permanent fact of life rather than a problem that appears now and then, optimization is less a specialized technique than a way of thinking that runs through nearly every applied field (Rajwar et al., 2023).Its reach is easy to underestimate. In engineering, optimization governs the sizing of structural members, the routing of transport and communication networks, and the calibration of manufacturing processes. In economics it underlies resource allocation, pricing policy, and the management of distribution systems. In artificial intelligence it drives model training, feature selection, and hyperparameter search; a clinical classifier for heart-disease prediction, for instance, gains measurable accuracy when its hyperparameters are tuned by a metaheuristic rather than fixed by hand(Hamajan, 2026). The same machinery reappears in sociology, in geology and resource exploration, and in genetics, where it supports sequence analysis and bioinformatic pipelines(Li et al., 2024). Everyday life differs only in stakes: deciding when to wake so that both sleep and deadlines are satisfied, or which route to take to work, are optimization decisions solved informally(Lang & Gao, 2025).What these cases share is a common structure, not a common subject. Once a problem has been written as a set of decision variables, a single measure of quality, and a set of constraints, the field it came from no longer matters. This is why optimization is sometimes called mathematical programming, and why a technique built for one field so often transfers unchanged to another. It is also why progress in optimization builds on itself: a single improvement in the underlying search machinery spreads at once to every field that uses the same formulation (Rajwar et al., 2023).

1.2. From Deterministic Methods to Metaheuristics

Optimization methods are usually divided into deterministic and heuristic families, with hybrid schemes drawing on both. A deterministic method follows a fixed sequence of steps and, given the same objective function and the same starting point, refollows the same path to the same answer every time. Heuristic methods add randomness on purpose: they give up the guarantee of identical results in exchange for the ability to climb out of the valley where the search happens to start (Darvishpoor et al., 2023).

For four decades the deterministic branch dominated engineering practice. Numerical linear and nonlinear programming — gradient descent, Newton and quasi-Newton schemes, sequential quadratic programming — are efficient, well understood, and on smooth convex models reliably reach the global optimum. The difficulty is that most real problems are none of those things. They need gradient information that may not exist; they contain many separate valleys, each with its own low point, so the answer returned depends on where the search began rather than on the problem itself; and where the objective or the constraints are rough and uneven, with sharp jumps, gradient-based procedures become unstable or fail clear (Tu et al., 2020). More detailed models have made this worse rather than better: as industrial and scientific models have grown in detail, obtaining a usable answer within an acceptable time has become harder, not easier (Li et al., 2024).

Metaheuristics were developed to work under exactly these conditions. They treat the objective as a black box (Boussaïd et al., 2013; Talbi, 2009a)— they only need to know the score a candidate solution receives, not how that score is computed. They require no derivatives, assume nothing about continuity or convexity, and combine fixed rules with random elements, often in a way borrowed from some natural or artificial process. Their appeal is as much practical as theoretical: they are simple to implement, easy to adapt to an unfamiliar problem, and ask almost nothing about how the model is built(Boussaïd et al., 2013).

Whatever metaphor a method adopts, its behavior is governed by two opposing forces. Exploration spreads the candidate solutions widely across the search space, searching for the region that holds the best overall answer; exploitation focuses effort around solutions already known to be good, refining them toward the best answer in their immediate neighborhood (Akbari et al., 2025). Randomness is what makes exploration possible, and it is also what allows a population to leave a local optimum once trapped in one. Too much exploration and search never stop moving long enough to converge; too much exploitation and it converges in the wrong place. Almost every design decision inside a metaheuristic is, in the end, a decision about how to move between these two modes, and when (Agrawal et al., 2021; Darvishpoor et al., 2023).

1.3. Research Activity in the Field

The result has been a rapid growth in numbers with few parallels in applied mathematics. More than five hundred distinct metaheuristics have now been published, a clear majority of them within the past decade, and the surrounding literature of surveys, comparative studies, and application papers has grown at a similar rate (Ameen, 2024; Rajwar et al., 2023; Shaikh et al., 2025). Research in the area falls into six recognizable streams (Ameen, 2024; Rajwar et al., 2023):

• Designing new algorithms, which accounts for the largest share of publications and for the taxonomy set out below. Five of the forty-two comparators illustrate the range. The Grey Wolf Optimizer guides a population using a ranking of leaders, with no velocity term (Mirjalili et al., 2014); Cuckoo Search uses Lévy-flight steps to widen the tail of the step-length distribution (Yang & Deb, 2009); the Honey Badger Algorithm alternates a digging phase with a honey-guided phase to switch between search modes (Hashim et al., 2022); the Archimedes Optimization Algorithm takes its operators from buoyancy and density rather than from animal behavior (Hashim et al., 2021); and the Circle Search Algorithm builds its update rule from the tangent relation on a circle, using no natural metaphor at all (Qais et al., 2022).
• Enhancement and modification of existing methods, where a known weakness is diagnosed and repaired. MRSO restructures the update rule of rat swarm optimization to correct an imbalance between its two search phases(Abdulla et al., 2024); EMRA revises the mud ring algorithm along comparable lines(Ali & Ameen, 2025); EMRSO adds hybrid memory strategies to MRSO (Hamid et al., 2026); EESB-FDO improves the fitness-dependent optimizer through a redesigned boundary-handling mechanism(Faraj et al., 2025); and DAOA replaces the static coefficients of the arithmetic optimization algorithm with dynamically varying ones(Khodadadi et al., 2022). The pattern is long established: L-SHADE adds linear population size reduction to SHADE and remains a standard baseline three decades after differential evolution itself(Tanabe & Fukunaga, 2014); a random-walk variant strengthens the exploration of the grey wolf optimizer (Gupta & Deep, 2019); an improved salp swarm algorithm sharpens its performance on feature selection (Hegazy et al., 2020); a chaotic dragonfly variant replaces random parameters with chaotic maps (Sayed et al., 2019); and an improved manta ray foraging optimizer was developed for multi-level image thresholding (Houssein et al., 2021).
• Hybridization, in which two methods are combined so that the strengths of one compensate for the weaknesses of the other. CDDO–HS, for example, couples child drawing development optimization with harmony search to improve information exchange within the population(Ameen, Rashid, & Askar, 2023). The pairing is usually chosen so that an exploratory method supplies breadth and an exploitative one supplies depth: PSO–GWO fuses the exploitation of particle swarm optimization with the exploration of the grey wolf optimizer (Şenel et al., 2019), and BGWOPSO carries the same pairing into binary feature-selection problems (Al-Tashi et al., 2019). Simulated annealing is a common partner for this reason, as in its hybrid with the whale optimization algorithm (Mafarja & Mirjalili, 2017). Earlier examples combine an ant colony with a genetic algorithm for protein function prediction (Nemati et al., 2009) and more recent ones pair a gradient-based optimizer with the grey wolf optimizer for human activity recognition (Helmi et al., 2021).
• Comparative analysis, which addresses the practical question of which algorithm to select for a given problem class. The scale of these studies has grown with the field; SFOA, for instance, was validated against one hundred competing optimizers(Zhong et al., 2025). The methodology of such comparisons has itself become a subject of study, with explicit guidelines now published on benchmark choice, baseline selection, parameter tuning, and statistical testing (LaTorre et al., 2021).
• Single- and multi-objective formulations, since real problems frequently involve several conflicting criteria and therefore have a set of best trade-off solutions rather than a single optimum. NSGA-II established the standard approach through elitist non-dominated sorting with a crowding-distance operator (Deb et al., 2002); MOEA/D instead decomposes the problem into scalar subproblems solved simultaneously (Q. Zhang & Li, 2007); MOPSO incorporates Pareto superiority into particle swarm optimization by way of an external repository (Coello et al., 2004); NSGA-III extends the framework to many objectives using a set of reference points (Deb & Jain, 2014); and the multi-objective grey wolf optimizer applies the same archive principle to a swarm method (Mirjalili et al., 2016)..
• Tutorials and surveys, which orient newcomers and periodically take stock of what has accumulated (Rajwar et al., 2023). Worked tutorials serve the same purpose at the level of a single algorithm, setting out its equations and iterations step by step so that a newcomer can reproduce it, as in the tutorial on child drawing development optimization (Ameen, Rashid, & Askar¹, 2023).

1.4. Ways of Classifying Metaheuristic Algorithms

The field has grown in several different directions at once, and no single scheme sorts it completely. Which classification is useful depends on the question being asked. A reader who wants to know where an idea came from needs a different scheme from one who wants to know how much memory the method uses, how many parameters must be tuned before it will run, or whether it can be given to a parallel machine. Eighteen axes recur in the literature, and they fall into four natural groups: where the idea came from, how the search is organized, how the method is written and tuned, and what it searches. Table 1 and Figure 1 summarizes all eighteen; the paragraphs below set out each group in turn.

• Where the idea came from. The oldest question is whether the designer borrowed the idea from the natural world. Methods that did divide into swarm intelligence, built on the collective behavior of many simple agents, and evolutionary computation, built on selection and inheritance. Methods that did not are grouped as non-nature-inspired, tabu search being the clearest case, since its forbidden-move list comes from an idea about memory rather than from any living thing (Blum & Roli, 2003). Within the nature-inspired group, a second axis asks which part of nature: authors divide this into three subgroups (Siddique & Adeli, 2015), four (Fister et al., 2013; Rajwar et al., 2023), five (Safiri & Nikoofard, 2023; Shaikh et al., 2025), or six (Khoshniat et al., 2024). A third axis, more recent and more critical, asks how faithfully a method follows the process it claims to copy. A strongly inspired algorithm genuinely reproduces the mechanism; a weakly inspired one uses the metaphor as decoration, and its equations could have been written without ever naming the animal in the title (Sörensen et al., 2025; Steer et al., 2009).
• How the search is organized. A second group of axes describes the machinery of the search itself. The most consequential counts how many candidate solutions are carried: one, as in tabu search, hill climbing and simulated annealing, or many, as in the genetic algorithm and ant colony methods (Jaddi & Abdullah, 2020). A second asks what the algorithm remembers. A memoryless method decides its next move from the current state alone and can look again at ground it has already covered; a method with memory keeps a record of earlier iterations and uses it to guide, which serves exploration and exploitation at the same time and usually pays for itself in fewer function evaluations. A third asks how a solution is assembled: an iterative method refines a candidate again and again, while a greedy method takes whatever looks best at each step and never revises it — faster and simpler, but it can easily miss the better answer that only becomes visible later (Sörensen et al., 2025). Two further axes, less often used but useful, ask whether the algorithm keeps one definition of which solutions count as neighbors or switches between several, as variable neighborhood search does, and whether it leaves the objective function alone or alters it during the run to push the search out of ground already visited, as guided local search does (Blum & Roli, 2003). A sixth asks how the next generation is actually produced: by shifting an existing solution, or by assembling a new one from several others, either through direct combination or through indirect coordination such as ant pheromone (Molina et al., 2020).

  

  Figure 1. Taxonomy of the 42 comparator algorithms by source of inspiration, after Rajwar et al. (2023): four families in the inner ring and one algorithm per segment in the outer ring, ordered by year within each family.

  Figure 1. Taxonomy of the 42 comparator algorithms by source of inspiration, after Rajwar et al. (2023): four families in the inner ring and one algorithm per segment in the outer ring, ordered by year within each family.

  

• How the method is written and tuned. A third group concerns the form of the algorithm rather than its behavior. Yang separates methods defined by a sequence of operations, such as evolution strategies and ant colony optimization, from those defined by explicit update equations, such as particle swarm optimization and cuckoo search; the great majority of modern methods are of the second kind, so what really separates one from another is the mathematical form of the update rule and not the set of terms of the metaphor (Yang, 2020). A newer axis simply counts the parameters a user must set before the method will run, grouping algorithms as free-parameter, mono-, bi-, tri-, tetra- and penta-parameter, with anything above five placed in a miscellaneous group; the fewer parameters a method carries, the less tuning it demands, which is why low-parameter methods are preferred for practical, real-world use (Rajwar et al., 2023). A related axis asks not how many parameters there are but how their values are set: on a fixed schedule, by feedback from the search, or by encoding them in the individuals so that they evolve with the population (Eiben et al., 1999).
• What is searched, and how it is deployed. A final group concerns the setting in which the method is used. A metaheuristic searches the space of solutions directly; a hyper-heuristic searches the space of heuristics, selecting or generating the low-level method that will then be applied, with generality rather than peak performance as the goal (Burke et al., 2013) .Methods also differ in whether they are defined as a single sequential process or as several processes running at once, since migration between subpopulations in an island model changes the search behavior and is not just a matter of how the code is written (Alba, 2005). They differ in the domain they natively address, continuous or binary, and a large share of published variants exist only to carry a continuous method into a binary setting by way of a transfer function (Pan et al., 2023). They differ in how they keep solutions within the rules. Some add a cost for breaking a rule; some use operators that cannot produce an illegal solution in the first place; some repair an illegal solution afterwards; and some treat rule-breaking as a separate measure to be reduced alongside the objective (Coello Coello, 2002). And they differ in how many criteria they optimize, a single-objective method returning one answer and a multi-objective method returns a set of trade-off solutions, called a Pareto front, in which no solution can be improved on one criterion without being made worse on another (Deb et al., 2002). Where two methods are combined, the combination itself can be classified by level, by order of execution, and by control strategy (Raidl, 2006; Talbi, 2009b).

  Table 1. Eighteen axes along which metaheuristic algorithms are classified. Every algorithm sits at a position on all eighteen simultaneously; the axes are complementary rather than competing.

  #	Axis	One pole	The other pole	Reference
  1	Source of the idea	nature-inspired	non-nature-inspired	(Blum & Roli, 2003)
  2	Domain of inspiration	three to six subgroups	—	(Rajwar et al., 2023)
  3	Strength of inspiration	strong	weak	(Steer et al., 2009)
  4	Candidates carried	single solution	population	(Jaddi & Abdullah, 2020)
  5	Use of memory	memory-based	memoryless	(Yasear & Ku-Mahamud, 2019)
  6	Assembly of a solution	iterative	greedy	(Sörensen et al., 2025)
  7	Neighborhood structures	one	several	(Blum & Roli, 2003)
  8	Objective function	static	dynamic	(Blum & Roli, 2003)
  9	Population movement	vector movement	solution creation	(Molina et al., 2020)
  10	Form of definition	procedure-based	equation-based	(Yang, 2020)
  11	Number of parameters	free parameter	up to penta- and beyond	(Rajwar et al., 2023)
  12	Parameter control	fixed schedule	adaptive or self-adaptive	(Eiben et al., 1999)
  13	What is searched	solutions	heuristics (hyper-heuristic)	(Burke et al., 2013)
  14	Execution	sequential	parallel or distributed	(Alba, 2005)
  15	Variable domain	continuous	binary or discrete	(Pan et al., 2023)
  16	Constraint handling	penalty	repair, feasible-only, separation	(Coello Coello, 2002)
  17	Number of objectives	single objective	multi-objective	(Deb et al., 2002)
  18	Hybrid structure	high-level, collaborative	low-level, integrative	(Raidl, 2006)

  Table 1. Eighteen axes along which metaheuristic algorithms are classified. Every algorithm sits at a position on all eighteen simultaneously; the axes are complementary rather than competing.

  #	Axis	One pole	The other pole	Reference
  1	Source of the idea	nature-inspired	non-nature-inspired	(Blum & Roli, 2003)
  2	Domain of inspiration	three to six subgroups	—	(Rajwar et al., 2023)
  3	Strength of inspiration	strong	weak	(Steer et al., 2009)
  4	Candidates carried	single solution	population	(Jaddi & Abdullah, 2020)
  5	Use of memory	memory-based	memoryless	(Yasear & Ku-Mahamud, 2019)
  6	Assembly of a solution	iterative	greedy	(Sörensen et al., 2025)
  7	Neighborhood structures	one	several	(Blum & Roli, 2003)
  8	Objective function	static	dynamic	(Blum & Roli, 2003)
  9	Population movement	vector movement	solution creation	(Molina et al., 2020)
  10	Form of definition	procedure-based	equation-based	(Yang, 2020)
  11	Number of parameters	free parameter	up to penta- and beyond	(Rajwar et al., 2023)
  12	Parameter control	fixed schedule	adaptive or self-adaptive	(Eiben et al., 1999)
  13	What is searched	solutions	heuristics (hyper-heuristic)	(Burke et al., 2013)
  14	Execution	sequential	parallel or distributed	(Alba, 2005)
  15	Variable domain	continuous	binary or discrete	(Pan et al., 2023)
  16	Constraint handling	penalty	repair, feasible-only, separation	(Coello Coello, 2002)
  17	Number of objectives	single objective	multi-objective	(Deb et al., 2002)
  18	Hybrid structure	high-level, collaborative	low-level, integrative	(Raidl, 2006)

These eighteen axes do not compete with one another. An algorithm has a position on every one of them at once: differential evolution is nature-inspired, evolutionary, population-based, memoryless, iterative, equation-based, bi-parameter, sequential by default and continuous by design. Three of the axes matter directly for what follows. The single- versus population-based distinction determines the class of method to which the present proposal belongs; the source of inspiration organizes the forty-two algorithms against which it is benchmarked; and the number of control parameters is the measure against which its own cost in tuning effort must later be assessed.

The single- versus population-based distinction is the most consequential for what follows. Single-solution methods such as simulated annealing and tabu search refine one candidate within its immediate neighborhood; they are economical and strongly exploitative, but a single trajectory can be captured by a local optimum with no mechanism for escape. Population-based methods instead maintain many candidates distributed across the search space, preserving diversity globally and sampling the landscape from several directions at once (Yang, 2020). This makes them markedly more robust on multimodal problems, and it is why essentially all recent proposals — including every algorithm considered in this study — are population-based (Jaddi & Abdullah, 2020).

1.5. Metaheuristics by Source of Inspiration

Within the population-based class, the most common organizing principle is the source of inspiration. It is also the oldest, and different authors partition it differently: into three subgroups (Siddique & Adeli, 2015), four (Fister et al., 2013), five (Safiri & Nikoofard, 2023; Shaikh et al., 2025), or six (Khoshniat et al., 2024; S. S. & H. S., 2022).

The present study adopts the four-subgroup scheme of Rajwar et al. (2023), which those authors identify as the most widely recognized of them, and which has the further value of dividing the physical-law family explicitly into a physics-based and a chemistry-based branch. The forty-two algorithms against which the present work is benchmarked are distributed across all four subgroups and are used below to illustrate each in turn; their distribution is shown in Figure 1, and their full details are given in Table 2.

• Evolutionary algorithms (EAs) are inspired by the Darwinian ideas of natural selection and the survival of those best suited to their surroundings. A population of individuals goes through simulated reproduction and mutation to create successive generations of offspring, the process being repeated so as to retain the genetic material that makes an individual better adapted to its environment while eliminating that which makes it weaker (Sloss & Gustafson, 2020). Darwin’s theory of natural evolution motivates both the genetic algorithm and differential evolution(Chakraborty et al., 2023), while genetic programming rests on the wider idea of biological evolution; gene expression programming, learning classifier systems, neuroevolutionary, and evolution strategies extend the same principle to other representations. Only one member of this family appears among the comparators. Differential Evolution remains the original model for the whole family and, three decades after its introduction, is still a serious competitor rather than a historical reference: its difference-vector mutation operator still sits at the heart of the best-performing entries in today’s optimization competitions (Tanabe & Fukunaga, 2014).
• Swarm intelligence (SI) algorithms model the collective behavior of systems that organize themselves with no central controller. Although Beni and Wang first used the term in 1989 in the context of cellular robotic systems, the term now covers any such system, and what marks these systems out is that they adapt as they go, communicate constantly, and share what they learn among agents that could neither defend themselves against a large hunting animal nor secure food alone (Rajwar et al., 2023). This is by a wide margin the largest family in literature, and it supplies twenty-seven of the forty-two comparators. Particle Swarm Optimization established the template, updating each particle from its own best position and the swarm’s (Abualigah et al., 2024; Kennedy & Eberhart, 1995); Cuckoo Search introduced Lévy-flight steps to widen the tail of the step-length distribution (Yang & Deb, 2009); and the Grey Wolf Optimizer later showed how a hierarchy of leaders could steer a population without an explicit velocity term (Mirjalili et al., 2014). Twenty-one further animal-behavior optimizers propose different mechanisms for regulating the same underlying trade-off: the Spotted Hyena Optimizer (Dhiman & Kumar, 2017), the Chimp Optimization Algorithm (Khishe & Mosavi, 2020), the Aquila Optimizer (Abualigah et al., 2021), the Honey Badger Algorithm (Hashim et al., 2022), the Snake Optimizer (Hashim & Hussien, 2022), the Mud Ring Algorithm (Desuky et al., 2022), the Zebra Optimization Algorithm (Trojovská et al., 2022), FOX (H. Mohammed & Rashid, 2023), the Gazelle Optimization Algorithm (Agushaka et al., 2023), Arctic Puffin Optimization (Wang et al., 2024), the Elk Herd Optimizer (Al-Betar et al., 2024), the Fossa Optimization Algorithm (Hamadneh et al., 2024), GOOSE (Hamad & Rashid, 2024), Modified Rat Swarm Optimization (Abdulla et al., 2024), the Starfish Optimization Algorithm (Zhong et al., 2025), the Scorpion Hunting Strategy (Singh et al., 2024), the Chinese Pangolin Optimizer (Guo et al., 2025), the Dhole Optimization Algorithm (B. O. Mohammed et al., 2025), the enhanced Mud Ring Algorithm (Ali & Ameen, 2025), Rüppell’s Fox Optimizer (Braik & Al-Hiary, 2025b), and Dwarf Sperm Whale Optimization (J. Zhang, 2026). A smaller non-swarm branch of the same bio-inspired family models plant rather than animal behavior: the Flower Pollination Algorithm (Yang, 2012), the Phototropic Growth Algorithm (Bohat et al., 2025), and an optimizer derived from water uptake and transport in plants (Braik & Al-Hiary, 2025a).
• Physical law-based algorithms (PhAs) encode physical or chemical law directly in the update rule, which often gives cleaner mathematics than a behavioral analogy. Rajwar (Rajwar et al., 2023) divide the family in two. Physics-based methods take gravitation, fields, the big bang, the black hole, and the galaxy as their primary sources, as in simulated annealing, the gravitational search algorithm, central force optimization, charged system optimization, electro-magnetism optimization, the big bang–big crunch algorithm, the black hole algorithm, and the galaxy-based search algorithm. Chemistry-based methods instead model molecular reaction, Brownian motion, and molecular radiation, as in gases Brownian motion optimization, the artificial chemical process, the ions motion optimization algorithm, and thermal exchange optimization. Four comparators fall in the physics-based branch: Transient Search Optimization, which borrows the response of a switched electrical circuit (Qais et al., 2020); the Archimedes Optimization Algorithm, whose operators derive from buoyancy and density (Hashim et al., 2021); the Young’s Double-Slit Experiment optimizer, which maps interference fringes onto search intensity (Abdel-Basset et al., 2023); and Mirage Search Optimization, which models atmospheric refraction (He et al., 2025). None of the forty-two takes up the chemistry-based branch. This family is of particular relevance here, since the present work draws on a physical system of the same kind.
• Miscellaneous is the residual subgroup, and in the formulation adopted here it collects methods founded on human behavior, game strategy, mathematical theorems, politics, artificial thought, and comparable ideas belonging to none of the preceding three (Rajwar et al., 2023). Ten comparators fall here, in two clearly separate groups. The first abstracts social, cognitive, or institutional processes. Harmony Search, modelled on a musician’s search for pleasing chords, is the long-standing representative (Zong Woo Geem et al., 2001); more recent entries include War Strategy Optimization (Ayyarao et al., 2022), the Skill Optimization Algorithm (Givi & Hubalovska, 2023), the Social Psychology Interaction Behavior Algorithm (Ameen, 2024), the Divine Religions Algorithm (Mozhdehi et al., 2025), and the Dream Optimization Algorithm (Lang & Gao, 2025). The Divine Religions Algorithm (Mozhdehi et al., 2025) is notable in the present context: it demonstrates that inspiration drawn from religious and social organization can yield a competitive optimizer and can be published in a well-regarded journal, provided the mathematics is stated independently of the metaphor. The second branch does away with natural inspiration altogether and constructs its operators from analytical or geometric objects. The Dynamic Arithmetic Optimization Algorithm builds its search on the four elementary arithmetic operators (Khodadadi et al., 2022); the Circle Search Algorithm exploits the tangent relation on a circle (Qais et al., 2022); the Weighted Average Algorithm updates from a weighted mean of the population (Cheng & De Waele, 2024); and the Newton Downhill Optimizer adapts a classical root-finding scheme to stochastic search (Xiao et al., 2026). Their advantage is transparency: the mechanism can be inspected directly, without decoding a metaphor first.

1.6. Limitations and Open Problems

Population-based methods are not without cost. They use far more memory and computation than single-solution alternatives, particularly at large population sizes or on expensive objective functions (Payerle et al., 2023). They fit naturally with parallel and distributed execution, which recovers part of that cost (Pardo et al., 2024), but they remain vulnerable to premature convergence: once population diversity collapses, the search stops making progress in whatever region it happens to be in, and most designs provide no mechanism for recovering the diversity that has been lost (Agrawal et al., 2021; Omeka et al., 2024).

These considerations set the requirements for any new proposal. It must balance exploration against exploitation rather than favoring either; it must maintain population diversity for long enough to avoid stagnation; and it must have reliably across problem classes rather than excelling on one family and failing badly on another. As well as these performance requirements, the field has begun to apply a further requirement about method. A substantial body of critical work argues that many recent metaheuristics contribute a new metaphor without contributing a new search mechanism — that once the terminology is removed, the update equations reduce to reparameterizations of methods already in the literature (Aranha et al., 2022; Dorigo, 2020; Sörensen & Glover, 2013).

The No Free Lunch theorem (Wolpert & Macready, 1997)supplies the formal justification for continuing to develop new methods at all. It establishes that, averaged over the space of all possible objective functions, no search algorithm outperforms any other: superior performance on one class of landscapes must be paid for by inferior performance on another (Velasco et al., 2024). The theorem therefore allows the search for algorithms whose inductive biases suit, practically important problem classes, while ruling out any claim of universal dominance. It is an argument for principled specialization, not for unlimited rapid growth in number (Velasco et al., 2024).

1.7. Motivation: Solar Geometry as a Control Schedule

Against that background, this work addresses a specific and, to our knowledge, unexamined weakness. Nearly every population-based metaheuristic requires some schedule for shifting emphasis from exploration to exploitation as the run proceeds, and in practice that schedule is almost always a hand-set value that only ever gets smaller — the linear decay of GWO, the quadratic and exponential decays used elsewhere, or a chaotic map substituted when additional variability is wanted. Such schedules share three defects. Being monotonic, they cannot re-expand the search radius once diversity has been lost. Their functional form and constants are chosen by the designer rather than derived, which adds tuning burden and obscures why one choice should be preferred to another. And they carry no physical meaning, so there is no principled basis for preferring any particular form.

Solar geometry offers an alternative that is free of all three defects, and it occupies a position that the taxonomy above leaves noticeably empty. Astronomical inspiration belongs formally to the physics-based branch of PhAs, whose recognized sources are gravitation, fields, the big bang, the black hole, and the galaxy; yet within that branch no established method draws on the observational geometry of the Earth–Sun system. Although optimizers have been derived from optics, electrical transients, buoyancy, and atmospheric refraction, and although cosmological objects have supplied metaphors for a generation of physics-based methods, the apparent motion of the Sun across the sky, and the shadow geometry that follows from it, have attracted almost no attention. For an observer at latitude $\phi$ on day $N$ of the year, the solar altitude angle $h\left(t\right)$ — obtained from the solar declination and hour angle by the standard spherical-astronomy relation — traces a smooth and strictly non-monotonic curve that rises from below the horizon, peaks at solar transit, and descends again. The associated shadow ratio, cot $h\left(t\right)$, behaves exactly as one would want a step-size coefficient to behave it diverges near sunrise and sunset, when the Sun is low and shadows are long, and contracts to a minimum at solar noon. A search radius driven by shadow length is therefore wide during the opening and closing situations of the cycle and narrow at its center, and it recovers automatically rather than decaying irreversibly. The whole schedule follows from $\phi$ and $N$ through closed-form trigonometric expressions containing no fitted constants, so an entire family of schedules is available at zero tuning cost.

The five canonical Islamic prayer times supply the natural partition of this curve. Fajr, Dhuhr, Asr, Maghrib, and Isha are not chosen by hand divisions of the day; each is fixed by an explicit geometric condition on the Sun’s position — Fajr and Isha by angle of the Sun below the horizon at twilight angles below the horizon, Dhuhr by the instant of solar transit, Asr by a prescribed ratio between an object’s shadow and its own length, and Maghrib by the disappearance of the solar disc corrected for atmospheric refraction. These conditions have been evaluated by standard astronomy for more than a thousand years, and they constitute a well-established segmentation of the solar cycle into five intervals of distinct geometric character: two low-altitude twilight regimes, one transit regime, one shadow-ratio regime, and one horizon-crossing boundary event. It is this partition, rather than the metaphor as such, that the proposed algorithm exploits.

1.8. Contributions and Organization of the Paper

The present study therefore introduces the Pray Algorithm Optimizer (PAO), a population-based metaheuristic whose control schedule and search radius are derived from solar altitude and shadow ratio, and whose search proceeds through five operators corresponding to the five prayer intervals. Across one solar day the controller switches on at most one extra operator per iteration, and the five together span the full range of search behaviors. Fajr restarts the worst individuals by reflecting them through the center of the search box to the point directly opposite, injecting long-range diversity that a contracted population can no longer produce on its own. Dhuhr draws a random subset part of the way toward the population median — a neighborhood statistic that settles the swarm without locking it to any single region. Asr and Maghrib both rely on global information, the best solution found so far, but to opposite ends: Asr takes small, distance-scaled Gaussian steps that refine individuals already close to the best, whereas Maghrib pulls its subset part of the way straight toward it, most strongly late in the run. Isha closes the cycle by giving stagnated individuals a heavy-tailed Lévy jump, restoring the diversity that the exploitative phases have used up. Underneath all five, the base engine of Section 3.4 supplies the two-at-a-time, difference-based mutation that runs on every iteration, so the phase operators add to rather than replace it.

The main contributions of this study are as follows:

• A new population-based metaheuristic, PAO, whose exploration–exploitation schedule and search radius are obtained in closed form from spherical-astronomy relations for solar altitude and shadow ratio, rather than from a hand-designed decay function. This places it in an inspiration category — astronomical — that is essentially unoccupied in the current literature.
• A non-monotonic, parameter-free control schedule that allows search diversity to be recovered in the later stages of a run, in contrast to the monotonically decreasing coefficients used by the majority of existing methods.
• Five structurally distinct search operators derived from the geometric conditions that define the five prayer intervals, together with an explicit term-by-term demonstration that the resulting update rules are not reparameterizations of PSO, DE, GWO, or SCA.
• A large-scale empirical evaluation against forty-two algorithms published between 1995 and 2026 and spanning all four inspiration families, on seventy-five benchmark functions: 23 classical functions, 30 CEC2017 functions, 10 CEC2019 functions, and 12 CEC2022 functions, under a uniform protocol of identical population size, identical evaluation budget, and 30 independent runs.
• A complete statistical treatment comprising 3,150 pairwise Wilcoxon rank-sum comparisons, Friedman rank analysis with post-hoc procedures, convergence and boxplot analysis, and a runtime comparison.
• A population-size sensitivity study over N ∈ {20, 30, 40, 50, 60} across all four benchmark suites.
• Application to 22 real-world engineering and application problems: Speed reducer (gearbox), Tension/compression spring, Pressure vessel, Three-bar truss, Gear train ratio, Cantilever beam, I-beam deflection, Tubular column, Piston lever, Car side impact, Welded beam, Reinforced concrete beam, Multi-disk clutch brake, Rolling element bearing, Belleville disc spring, Hydrostatic thrust bearing, Robot gripper, Gas transmission compressor, Stepped cone pulley, Planetary gear train, 10-bar planar truss (FEM), Heat exchanger (Himmelblau)
• Public release of the PAO source code at [ https://www.mathworks.com/matlabcentral/fileexchange/184442-pao-for-global-and-constrained-engineering-optimization ], so that the results reported here can be reproduced independently.

The remainder of this paper is organized as follows. Section 2 develops the astronomical background, including the computation of solar position and of the five prayer times. Section 3 presents PAO, its mathematical model, pseudocode, complexity analysis, and its relationship to existing algorithms. Section 4 describes the experimental protocol. Section 5 reports and discusses the benchmark results. Section 6 applies PAO to real-world engineering problems. Section 7 discusses limitations, and Section 8 concludes.

2. Astronomical Background and Inspiration

This section explains the idea behind the algorithm before any mathematics is applied to it. It first describes the five daily prayers and the fact that their times are not chosen by people but fixed by the position of the Sun. It then shows how the Sun’s position can be written as a single angle that moves through one full turn each day, and how each prayer time corresponds to a particular value of that angle. The key point for the reader is that this gives a natural, repeating schedule that changes with the season and with where you are on Earth — exactly the kind of schedule an optimizer needs when it has to switch between wide searching and careful refining.

2.1. Prayer in Islam and the Determination of Prayer Times

Prayer (ṣalāh) is one of the five main practices of Islam, and the five daily prayers are something every adult Muslim is required to perform. Like the other main religious duties — paying zakat and giving voluntary charity (ṣadaqah), and dividing the property a person leaves behind after death — prayer is not left to each person to decide. Its number, the way it is performed, and above all its times are set by Islamic law (the Sharia) through its two main sources, the Qur’an and the Sunnah (Bahri & Hasibuan, 2024; Esposito, 2011).

That the five daily prayers must be performed at fixed times is a rule from the holy texts, not just a custom: the Qur’an describes prayer as a duty that must be done at set times, and the Prophetic tradition (Sunnah) gives, for each of the five, a starting point marked by a visible position of the Sun. So the rules for when each prayer may be performed belong to Islamic religious law (fiqh) just as much as to astronomy, and the two meet exactly at the point where a religiously defined moment has to be turned into a clock time (Rojak & Fawzi, 2024).

Each prayer begins when the Sun reaches a stated position — sitting a fixed distance below the horizon before sunrise, crossing the highest point of the sky, shadows growing to a set length, the Sun sinking below the horizon, and the evening light fading after sunset. Because these are astronomical events, the prayer times change from place to place and from day to day. Three things move them: where the observer is on Earth, given by latitude and longitude; the date, and so the time of year (season), through the Sun’s changing declination; and the rules chosen by the different schools of Islamic law (Bahri & Hasibuan, 2024; Meeus, 1991; Rojak & Fawzi, 2024).

The four Sunni schools of law — Hanafi, Shafi’i, Maliki, and Hanbali — agree on the astronomical events that mark the five prayers but differ on some of the exact points used to read them. The clearest difference is the afternoon prayer, ʿAṣr: the Shafi’i, Maliki, and Hanbali schools start it when an upright object’s shadow has grown by an amount equal to the object’s own height beyond its noon length, while the Hanafi school waits until the added shadow equals twice that height. There are also smaller differences in how far below the horizon the Sun must be at dawn and dusk (the twilight angles) for Fajr and Isha. These differences are why published prayer timetables for the same city on the same day do not always match, and why any calculation must state which rule it follows (Rojak & Fawzi, 2024).

What matters for this work is that the five prayer moments split the whole day using geometry alone: once we know a latitude, a date, and a chosen school’s rule, their positions come straight from the position of the Sun, with no free number left to adjust by hand. This is exactly the property PAO reuses — a smooth, repeating, up-and-down schedule set entirely by astronomy — turning it into a control schedule for optimization, as the rest of this section explains.

2.2. The Idea in Plain Terms

Every population-based optimizer needs a rule that decides, at each moment of a run, how far its candidate solutions should move. Early on, long moves are wanted, because the search must cover ground; later, short moves are wanted, because the best solution found so far needs fine-tuning rather than replacing. In almost all existing algorithms this rule is a hand-made number that only ever decreases.

PAO takes the rule from somewhere else: the position of the Sun during a single day. The Sun rises, climbs to its highest point, and sets. The shadow it casts is long at the two ends of the day and short in the middle. That pattern has four properties that are useful here. It is smooth. It is not monotone, so the quantity driven by it can grow again after getting smaller. It repeats exactly. And it is fixed completely by two numbers — the latitude of the observer and the day of the year — with no constant left for the designer to choose.

The five daily Islamic prayers divide this curve at instants that are themselves defined geometrically rather than by convention: two are fixed by how far the Sun sits below the horizon, one by the instant the Sun crosses the meridian, one by the ratio of an object’s shadow to its own length, and one by the Sun touching the horizon. The result is a partition of the day into intervals of distinct geometric character. PAO gives each interval its own search operator.

2.3. Where the Sun is

Let $n$be the day of the year $(1 to 365)$ and let $\phi$ be the latitude of the observer, positive north. The declination $\delta$ is the angle between the Sun’s rays and the plane of the Earth’s equator; it varies slowly through the year and is given to sufficient accuracy by Cooper’s equation (Cooper, 1969; Duffie & Beckman, 2013):

$$\mathit{\delta }={\mathbf{23.45}}^{°}\times \mathit{sin}\left(\left({\displaystyle \frac{\mathbf{360}}{\mathbf{365}}}\right)\times \left(\mathbf{284}+\mathit{n}\right)\right)$$

(1)

Local time is measured by the hour angle $\omega$, which is zero at solar noon and advances by 15° for every hour. Writing $t$ for clock time in minutes after midnight (Duffie & Beckman, 2013; Ma & Iqbal, 1983),

$$\mathit{\omega }\left(\mathit{t}\right)={\mathbf{15}}^{°}\times \left({\displaystyle \frac{\mathit{t}}{\mathbf{60}}}-\mathbf{12}\right)$$

(2)

The height of the Sun above the horizon, the solar altitude $\alpha$, then follows from a standard astronomy formula(Duffie & Beckman, 2013; Meeus, 1991):

$$\mathbf{sin}\mathit{\alpha }=\mathbf{sin}\mathit{\phi }\times \mathbf{sin}\mathit{\delta }+\mathbf{cos}\mathit{\phi }\times \mathbf{cos}\mathit{\delta }\times \mathbf{cos}\mathit{\omega }$$

(3)

Two quantities derived from $\alpha$ are used later. The zenith angle measures the Sun’s distance from straight overhead, and the shadow ratio is the length of the shadow cast by a vertical object of unit height (Ma & Iqbal, 1983):

$$\mathit{z}={\mathbf{90}}^{°}-\mathit{a}$$

(4)

$$\mathit{s}\left(\mathit{t}\right)=\mathbf{cot}\mathit{a}\left(\mathit{t}\right)$$

(5)

The shadow ratio is the quantity that behaves like a search radius. It becomes very large when the Sun is low, near sunrise and sunset, and reaches its minimum at solar noon. A step-size driven by it is therefore wide at the two ends of the cycle and narrow in the middle, and it widens again on its own rather than decaying once and for all.

2.4. Reading the Geometry Backwards

Prayer times are defined by the Sun reaching a stated height, so $Eq. \left(3\right)$ has to be solved for $\omega$ rather than evaluated. Rearranging gives the hour angle at which the altitude equals a set prescribed value α (Meeus, 1991):

$$\mathbf{cos}\mathit{\omega }\left(\mathit{\alpha }\right)={\displaystyle \frac{\mathbf{sin}\mathit{\alpha }-\mathbf{sin}\mathit{\phi }\times \mathbf{sin}\mathit{\delta }}{\mathbf{cos}\mathit{\phi }\times \mathbf{cos}\mathit{\delta }}}$$

(6)

Sunrise and sunset are the case $\alpha = -0.833°$, the small negative value accounting for atmospheric refraction and the radius of the solar disc. With $\omega s = \omega (-0.833°)$ the two instants and the length of the day are (Meeus, 1991):

$$\mathit{D}\mathit{d}\mathit{a}\mathit{y}=\mathit{t}\mathit{s}\mathit{s}-\mathit{t}\mathit{s}\mathit{r}, \mathit{w}\mathit{h}\mathit{e}\mathit{r}\mathit{e}\left\{\begin{array}{c}\mathbf{t}\mathbf{s}\mathbf{r}=\mathbf{12}-{\displaystyle \frac{\mathit{\omega }\mathit{s}}{\mathbf{15}}}\\ \mathbf{t}\mathbf{s}\mathbf{s}=\mathbf{12}+{\displaystyle \frac{\mathit{\omega }\mathit{s}}{\mathbf{15}}}\end{array}\right.$$

(7)

2.5. The Five Prayer Times as Geometric Conditions

Each prayer begins when the Sun satisfies one specific geometric condition. Table 3 states them. Only Asr needs an extra formula: its condition is that the shadow of a vertical object has grown by $k$ times the object’s own height beyond its length at noon, with $k = 1$ in the Shafi‘i, Maliki and Hanbali schools and $k = 2$ in the Hanafi school. Converting that statement into an altitude gives(Bahri & Hasibuan, 2024; Rojak & Fawzi, 2024):

$$\mathbf{cot}\mathit{\alpha Asr}=\mathit{k}+\mathbf{tan}|\mathit{\phi }-\mathit{\delta }|$$

(8)

so $\alpha Asr$ follows immediately, and the corresponding instant follows from $Eq. \left(6\right)$. In general, every prayer instant is obtained by substituting its altitude condition into $Eq. \left(6\right)$:

$$\mathit{t}\mathit{p}\mathit{r}\mathit{a}\mathit{y}\mathit{e}\mathit{r}=\mathbf{12}\pm {\displaystyle \frac{\mathit{\omega }\left(\mathit{\alpha }\mathit{c}\mathit{o}\mathit{n}\mathit{d}\mathit{i}\mathit{t}\mathit{i}\mathit{o}\mathit{n}\right)}{\mathbf{15}}}$$

(9)

with the minus sign before solar noon and the plus sign after it. Here $\alpha$_condition is not a fixed number but a placeholder for the solar altitude that defines the prayer in question: the value in the “Altitude used” column of Table 3 — $-18°$ for Fajr, $-17°$for Isha, and $-0.833°$ for Maghrib — and, for Asr, the altitude from $Eq. \left(8\right)$. Dhuhr is the exception: it is fixed by the Sun crossing the meridian ($\omega = 0$), not by a depression angle, so it is read directly rather than through $Eq. \left(9\right)$.

2.6. A Clock that Lengthens with the Season

The five instants of Table 3 move as $\phi$ and $n$ change: near the equator they shift little, at high latitude they shift a great deal. If the algorithm referred to them by clock time, its schedule would have to be re-derived for every location and date. To avoid this, time is re-expressed as an angle that always completes one full turn over one day (see Figure 2), with sunrise at $0°$ and sunset at $180°$. This remapping is introduced in the present work (Eqs. 8–9 are not a standard astronomical relation). Call it the day angle $\theta$:

$$\mathit{\theta }\left(\mathit{t}\right)={\displaystyle \frac{{\mathbf{180}}^{°}\times \left(\mathit{t}-\mathit{t}\mathit{s}\mathit{r}\right)}{\mathbf{1440}-\mathit{D}\mathit{d}\mathit{a}\mathit{y}}},\mathit{f}\mathit{o}\mathit{r} \mathit{d}\mathit{a}\mathit{y}\mathit{l}\mathit{i}\mathit{g}\mathit{h}\mathit{t},\mathit{t}\mathit{s}\mathit{r}\le \mathit{t}\le \mathit{t}\mathit{s}\mathit{s}$$

(10)

$$\mathit{\theta }\left(\mathit{t}\right)={\displaystyle \frac{{{\mathbf{180}}^{°}+\mathbf{180}}^{°}\times \left(\mathit{t}-\mathit{t}\mathit{s}\mathit{s}\right)}{\mathbf{1440}-\mathit{D}\mathit{d}\mathit{a}\mathit{y}}},\mathit{f}\mathit{o}\mathit{r} \mathit{n}\mathit{i}\mathit{g}\mathit{h}\mathit{t}$$

(11)

The inverse map, which converts an angle back to a clock time, is used to lay the five windows onto the run. Within daylight, $\theta$ is an affine function of the hour angle $\omega$ of $Eq. \left(2\right)$, so it carries the same information; the piecewise definition simply forces night to occupy the same $180°$ as day, whatever the season. The effect is that a window written in $\theta$ automatically stretches in summer and contracts in winter, and adapts to latitude, without a single number being changed.

Figure 3. The day-angle representation θ that PAO uses to place the five prayer windows. Daylight is mapped onto θ ∈ [0°, 180°] (θ = 0° at sunrise on the east, θ = 180° at sunset on the west) and night onto θ ∈ [180°, 360°], so a single angle labels every instant of the day. The windows are then read off as arcs of this circle: the forenoon gap up to θ ≈ 95° carries no prayer, followed by Dhuhr, Asr near θ ≈ 135°, Maghrib at θ = 180°, Isha from θ ≈ 200°, and Fajr near θ ≈ 340° — the same boundaries collected in the window table of section2.7. Because θ comes from Eqs. (10)–(11) rather than from clock time, these edges keep their meaning at any latitude or date.

Figure 3. The day-angle representation θ that PAO uses to place the five prayer windows. Daylight is mapped onto θ ∈ [0°, 180°] (θ = 0° at sunrise on the east, θ = 180° at sunset on the west) and night onto θ ∈ [180°, 360°], so a single angle labels every instant of the day. The windows are then read off as arcs of this circle: the forenoon gap up to θ ≈ 95° carries no prayer, followed by Dhuhr, Asr near θ ≈ 135°, Maghrib at θ = 180°, Isha from θ ≈ 200°, and Fajr near θ ≈ 340° — the same boundaries collected in the window table of section2.7. Because θ comes from Eqs. (10)–(11) rather than from clock time, these edges keep their meaning at any latitude or date.

2.7. The Five Windows, and how Close They are to the Real Thing

PAO fixes the five phases as intervals of $\theta$, given in Table 4. One interval of the day carries no phase: from sunrise to Dhuhr, $\theta \in (0°, 95°)$. This is done on purpose, and it matches the source of the metaphor, since no obligatory prayer falls in the forenoon. During that stretch the base engine of Section 3.4 runs unmodified.

Because $\theta$ is only piecewise linear in time, the window boundaries do not fall exactly on the true prayer instants. Table 5 quantifies the difference for a mid-latitude spring day. The largest error is 19 minutes out of a 24-hour cycle, that is 1.3% of the period, which is well inside the spread between the calculation conventions of the different juristic schools. The approximation is therefore adequate for the purpose it serves.

2.8. What the Optimizer Gains from this

The astronomy of Section 2.1, Section 2.2, Section 2.3, Section 2.4, Section 2.5, Section 2.6 and Section 2.7 turns out to supply exactly what a search schedule needs: a cycle that repeats, windows of unequal length, and an order that is not monotone. Each of these has a direct algorithmic counterpart — the day angle drives the iteration clock, the window boundaries select the active operator, and the cycle weight sets its strength. Table 6 states the correspondence term by term.

Two consequences are worth stating explicitly, because they are what distinguishes this schedule from a conventional decay. First, the schedule is not monotone: an exploratory phase (Isha) occupies the final third of every run, so population diversity can be recovered late rather than only lost. Second, the schedule contains no fitted constant. Changing $\phi$ or $n$ produces a different but equally valid schedule, which makes the family of schedules available for free and turns a design choice into an experimental variable (Figure 3).

Figure 5. The five phase operators and the solar geometry that generates them, for φ = 35°N on day 105. (a) Solar altitude over 24 hours, with the five prayer windows shaded and the twilight thresholds that define Fajr and Isha marked; the black triangles are the exact prayer instants of Table 3. (b) The day angle θ of Eq. (10)–(11), which maps daylight onto [0°, 180°] and night onto [180°, 360°]. (c) The control parameters that the active phase applies to the base engine. Generated by Figure 2_pao_phases.m; passing a different latitude or day of year regenerates the whole schedule.

Figure 5. The five phase operators and the solar geometry that generates them, for φ = 35°N on day 105. (a) Solar altitude over 24 hours, with the five prayer windows shaded and the twilight thresholds that define Fajr and Isha marked; the black triangles are the exact prayer instants of Table 3. (b) The day angle θ of Eq. (10)–(11), which maps daylight onto [0°, 180°] and night onto [180°, 360°]. (c) The control parameters that the active phase applies to the base engine. Generated by Figure 2_pao_phases.m; passing a different latitude or day of year regenerates the whole schedule.

2.9. Modelling Assumptions and the Reference Configuration

The astronomical schedule described above is exact in principle, but the reference implementation realizes it under three simplifications made on purpose. They are stated here so that the results reported later can be reproduced and the scope of the model understood.

First, the Earth’s surface at the observer is idealized as a flat, unobstructed plane. The observer is taken to stand on a level horizon, and no correction is applied for elevation above sea level or for local topography such as the surrounding mountain ranges; the dip of the horizon that arises from the Earth’s curvature and the observer’s height is likewise neglected. This is the standard idealization used in elementary prayer-time computation. It leaves the celestial geometry of Eqs. (1)–(6) unchanged — the Sun’s altitude is still obtained from the standard astronomy formula — and affects only the reference level of the local horizon against which sunrise, sunset, and the twilight angles are measured. In adopting it we set aside the actual elevation of the reference city (Chamchamal lies a few hundred meters above sea level, in a plain flanked by mountains); the resulting shift in the computed instants is small relative to the spread already introduced by the choice of juristic convention, and it does not alter the phase structure that PAO extracts from the schedule (Bahri & Hasibuan, 2024; Meeus, 1991).

Second, a single reference location is fixed for all worked examples and for the default clock of the implementation: the city of Chamchamal, in Sulaymaniyah Governorate, Kurdistan Region, Iraq, at approximately $35.5° N, 44.8° E$. This latitude is deliberately close to the $\phi = 35° N$ used in the worked example of Section 2.7, so that the validation figures reported there correspond to the reference city itself. Because the schedule takes latitude and day of year as its only inputs, however, nothing in the method depends on this choice: substituting the coordinates of any other location, or any other date, regenerates a complete and equally valid schedule at no additional tuning cost (see Figure 4).

Third, among the four Sunni schools the Shafi’i school of jurisprudence is adopted throughout. In practice this fixes a single quantity, the shadow factor of the afternoon prayer: the Shafi’i criterion corresponds to $k = 1$ in Eq. (10), so ʿAṣr begins when the shadow of an upright object exceeds its noon length by one object-height. The Hanafi value k = 2 would move only the ʿAṣr boundary and would leave the four remaining phases, and the overall exploration–exploitation profile of the schedule, unchanged. For Fajr and Isha the twilight depressions of $-18°$ and $-17°$ stated in Table 3 are used, consistent with the conventions commonly applied at this latitude (Rojak & Fawzi, 2024).

Figure 7. The prayer schedule projected onto the whole Earth at a single instant, for the Chamchamal reference configuration (spring, day 80, local 15:45). (a) Equirectangular map and (b) orthographic globe centered on Chamchamal; each color marks the region of the Earth in which a given prayer is in force at that moment, with the sub-solar point (yellow disc) and Chamchamal (triangle) marked. The figure makes concrete why the schedule must depend on latitude, longitude and date: at one instant different places sit in different prayers, and the local clock times listed in the inset (Fajr 04:33, Sunrise 05:58, Dhuhr 12:01, Asr 15:27, Maghrib 18:04, Isha 19:24) follow directly from the geometry of Section 2 rather than from a stored timetable.

Figure 7. The prayer schedule projected onto the whole Earth at a single instant, for the Chamchamal reference configuration (spring, day 80, local 15:45). (a) Equirectangular map and (b) orthographic globe centered on Chamchamal; each color marks the region of the Earth in which a given prayer is in force at that moment, with the sub-solar point (yellow disc) and Chamchamal (triangle) marked. The figure makes concrete why the schedule must depend on latitude, longitude and date: at one instant different places sit in different prayers, and the local clock times listed in the inset (Fajr 04:33, Sunrise 05:58, Dhuhr 12:01, Asr 15:27, Maghrib 18:04, Isha 19:24) follow directly from the geometry of Section 2 rather than from a stored timetable.

3. The Proposed PAO Algorithm

This section turns the idea of Section 2 into a working method. It maps a complete run of the algorithm onto a single day, so that each iteration falls at a particular clock time and therefore inside one of the five prayer windows. Each window switches on its own search behavior: Fajr restarts the worst solutions, Dhuhr pulls the group together, Asr makes small careful improvements, Maghrib moves strongly toward the best solution found so far, and Isha gives stuck solutions a large jump. The section gives the equations for each of these steps, the settings they use, the pseudocode and flowchart of the complete method, and an estimate of how much computing work it needs.

3.1. Problem Statement and Notation

PAO addresses the bound-constrained single-objective minimization problem:

$$\mathit{m}\mathit{i}\mathit{n}\mathit{i}\mathit{m}\mathit{i}\mathit{z}\mathit{e} \mathit{f}\left(\mathit{x}\right),\mathit{x}=\left({\mathit{x}}_1,...,{\mathit{x}}_{\mathit{d}}\right)\in {\mathcal{R}}^{\mathit{d}}, \mathit{l}\mathit{b}\mathit{j}\le \mathit{x}\mathit{j}\le \mathit{u}\mathit{b}\mathit{j}$$

(12)

Equation (12) states the optimization problem that PAO addresses: to locate the vector x that minimizes the objective function $f$. Here x = (x₁, …, x_d) is a candidate solution composed of d decision variables, ℝ^d denotes the space of real-valued vectors of length $d$, and $lb\_j$ and $ub\_j$ are the lower and upper bounds on the $j-th$ variable, so that each coordinate remains within its permitted interval. This is the standard formulation of a bound-constrained minimization problem and is adopted by PAO without modification. Where $f$is treated as a black box: no derivative, no continuity and no convexity is assumed. Table 7 lists the symbols used throughout.

3.2. How PAO is Put Together

PAO has two layers, and it is important to be clear about which is which:

The lower layer is the base engine: an adaptive, population-based search engine that mutates each candidate using a scaled difference between population members, keeps a small external archive of replaced solutions, and self-tunes its own step size and mixing rate as the run proceeds. It is a well-known, off-the-shelf method that we do not claim as new; it runs on every iteration and does most of the search.

The upper layer is the astronomical controller, and this is what PAO contributes. On each iteration it converts the iteration counter into a position in the solar cycle, decides which of the five prayer windows that position falls in, and then does two things: it multiplies the engine’s two control parameters by phase-specific factors, and it switches on at most one additional operator chosen to suit the character of that window. When the position falls in the forenoon gap, no additional operator runs at all.

3.3. Initialization

The starting solutions are scattered evenly, at random, across the whole search region — each coordinate is set to $x = lb + r\times (ub - lb)$ with $r$ drawn uniformly from $[0, 1]$ — and the best individual is recorded.

The two success-memories — the small tables that remember which step-size and mixing settings have worked well recently — are set to $MF = MCR = 0.5$, the archive — the running store of parents that have just been replaced — starts empty, and the counter that records how long each individual has gone without improving is reset to zero. The daily prayer-time schedule of Section 2 is worked out once before the loop, because it depends only on the latitude $\phi$ and the day of the year $n$.

3.4. The Base Engine

The lower layer is the adaptive, population-based base engine described in Section 3.2. Because it is established machinery and not part of the contribution, only it’s one characteristic step is written out; the remaining textbook steps are described in words. On each iteration, for every individual, the engine samples a scale factor F and a crossover rate $CR$ from two small memories of recently successful values — $F$ from a heavy-tailed random draw (a Cauchy distribution), $CR$ from a normal, bell-curve draw — and then shifts them by the phase multipliers $E_{scale}$ and ${CR}_{bias}$ supplied by the controller of Section 3.5, which is the only point at which the controller reaches into the engine. The new trial solution combines the current solution with a move toward one of the best $pN$ individuals and a scaled difference between a population member and a partner drawn from the population together with the archive:

$${\mathit{v}}_{\mathit{i}}={\mathit{x}}_{\mathit{i}}+{\mathit{F}}_{\mathit{i}}\times \left[\left({\mathit{x}\mathit{p}}_{\mathit{b}\mathit{e}\mathit{s}\mathit{t}}-{\mathit{x}}_{\mathit{i}}\right)+\left({\mathit{x}\mathit{r}}_{\mathbf{1}}-{\widetilde{\mathit{x}\mathit{r}}}_{\mathbf{2}}\right)\right]$$

(13)

Here ${xp}_{best}$is one of the current best individuals, selected at random from the top-ranked fraction; ${xr}_1$is a randomly chosen population member; ${\widetilde{xr}}_2$is a member drawn at random from the union of the population and the external archive of recently replaced solutions; and $F_i$ is the individual scaling factor. The base engine of Section 3.4 supplies this operator; PAO leaves its form unchanged and only rescales $F_i$according to the active phase.

The trial is then mixed with the parent coordinate by coordinate, always keeping at least one coordinate from the trial so it can never be an exact copy; the trial is brought back inside the box and replaces the parent only when it is at least as good. Every replaced parent is kept in an external archive of capacity $N$, and at the end of each iteration the two memories are refreshed from the $F$ and $CR$ values that produced improvements, each weighted by how much it improved the objective ((a special average that leans toward the larger successful values, for $F$, a weighted average for $CR$). These are all standard steps of the base engine and are not written out here.

The slot index $k$ advances by one each iteration and wraps around at $H$, so the memories hold a moving record of what has recently worked.

3.5. The Astronomical Controller

The run is mapped onto one solar day. Iteration it is placed at:

$$\mathit{t}=⌊{\displaystyle \frac{\mathbf{1440}\times \left(\mathit{i}\mathit{t}-\mathbf{1}\right)}{{Max}_{Iter}}}⌋\mathit{m}\mathit{i}\mathit{n}\mathit{u}\mathit{t}\mathit{e}\mathit{s} \mathit{a}\mathit{f}\mathit{t}\mathit{e}\mathit{r} \mathit{m}\mathit{i}\mathit{d}\mathit{n}\mathit{i}\mathit{g}\mathit{h}\mathit{t}$$

(14)

Equation (14) maps the iteration index onto a time of day, so that the prayer schedule can govern the search. Here it is the current iteration, ${Max}_{Iter}$ is the total number of iterations, 1440 is the number of minutes in a day, and ⌊·⌋ denotes the floor (round-down) operation. Consequently, the first iteration corresponds to 00:00 and the final iteration to approximately 24:00. This mapping is introduced in the present work and is the mechanism by which the astronomical schedule is connected to the optimizer.

The day angle $\theta$ follows from $Eq. \left(10\right)-\left(11\right)$. Comparing $\theta$ with the windows of Table 4 selects the phase. A weight that falls linearly from one to zero across the cycle, measured from the start of Fajr, controls how strongly the phase operators act:

$${\mathit{w}}_{\mathit{t}}=\mathbf{1}-{\displaystyle \frac{\left(\mathit{t}-{\mathit{t}}_{\mathit{F}\mathit{a}\mathit{j}\mathit{r}}\right)\mathit{m}\mathit{o}\mathit{d}\mathbf{1440}}{\mathbf{1440}}},{\mathit{w}}_{\mathit{t}}\in \left[\mathbf{0,1}\right]$$

(15)

Equation (15) defines the cycle weight $w$, a quantity that decreases from one to zero and indicates how far the run has advanced through the day. Here $t$ is the time obtained from Equation (14), ${\mathit{t}}_{\mathit{F}\mathit{a}\mathit{j}\mathit{r}}$is the start time of Fajr, the modulo-1440 operation wraps the elapsed time around the 24-hour clock, and division by 1440 scales the result to the interval [0, 1]. The weight therefore equals one at Fajr and decreases linearly to zero across the cycle. It is a design choice for the present work: the schedule is read from the solar day rather than from a fixed decay curve.

Operators that should be strong early use $w$ directly; the one that should strengthen late uses $1 - w$. Table 8 collects the phase settings.

3.6. Fajr — Opposition-Based Re-Seeding

Fajr occupies the darkest part of the cycle, just before first light, and its role is exploration. The worst $10\%$ of the population are reflected through the center of the search box to the point directly opposite, and the reflection is kept only if it is better:

$${\tilde{\mathit{x}}}_{\mathit{i}}=\mathit{l}\mathit{b}+\mathit{u}\mathit{b}-{\mathit{x}}_{\mathit{i}}, {\mathit{x}}_{\mathit{i}}=\left\{\begin{array}{c}{\tilde{\mathit{x}}}_{\mathit{i}}, if f\left({\tilde{\mathit{x}}}_{\mathit{i}}\right)<f\left({\mathit{x}}_{\mathit{i}}\right)\\ {\mathit{x}}_{\mathit{i},\mathit{o}\mathit{t}\mathit{h}\mathit{e}\mathit{r}\mathit{w}\mathit{i}\mathit{s}\mathit{e}}\end{array}\right.$$

(16)

Here ${\tilde{\mathit{x}}}_{\mathit{i}}=\mathit{l}\mathit{b}+\mathit{u}\mathit{b}-{\mathit{x}}_{\mathit{i}}$ is the opposite (mirror) point of ${\mathit{x}}_{\mathit{i}}$, and the reflected point is retained only if it improves the objective, $\mathit{f}\left({\tilde{\mathit{x}}}_{\mathit{i}}\right)<\mathit{f}\left({\mathit{x}}_{\mathit{i}}\right)$; otherwise ${\mathit{x}}_{\mathit{i}}$is left unchanged.

This is opposition-based learning restricted to the individuals that have least to lose. It costs one evaluation per re-seeded individual and can move a solution across the entire domain in a single step, which no differential mutation can do once the population has contracted.

3.7. Dhuhr — Median Pull

Dhuhr begins just after the Sun crosses the meridian, when shadows are shortest. Its role is contraction, but toward the population itself rather than toward the current best. Let $m$ be the coordinate-wise median of the population. A random 30% of individuals are drawn part of the way toward it:

$${\mathit{x}}_{\mathit{i}}\leftarrow {\mathit{x}}_{\mathit{i}}+{\mathit{\eta }}_{\mathit{D}}\times \mathit{r}\times \left(\mathit{m}-{\mathit{x}}_{\mathit{i}}\right),\mathit{r}~\mathit{U}\left(\mathbf{0,1}\right)$$

(17)

Here m is the coordinate-wise median of the population (a central value that is unaffected by extreme points), $r$ is a uniform random number in [0, 1], and ${\mathit{\eta }}_{\mathit{D}}$ = $0.15 w$ controls the step length, which diminishes as the day advances.

Moving toward the middle of the group rather than toward the best solution is what makes this a balancing step instead of a second exploitation step. The middle value is not affected much by a few unusual points, so the operator draws the population closer together without locking it onto whichever region currently holds the best solution.

3.8. Asr — Local Gaussian Refinement

Asr is defined by shadows lengthening past a fixed multiple of the object height, and it is the exploitation phase. A random 35% of individuals take a Gaussian step whose width in each coordinate is proportional to how far that coordinate currently sits from the best solution:

$${\mathit{z}}_{\mathit{i}}={\mathit{x}}_{\mathit{i}}+{\mathit{\eta }}_{\mathit{A}}\times \mathcal{N}\left(\mathbf{0,1}\right)\odot \left|{\mathit{x}}^{\mathit{*}}-{\mathit{x}}_{\mathit{i}}\right|,{\mathit{x}}_{\mathit{i}}=\left\{\begin{array}{c}{\mathit{z}}_{\mathit{i}}, if f\left({\mathit{z}}_{\mathit{i}}\right)<f\left({\mathit{x}}_{\mathit{i}}\right)\\ {\mathit{x}}_{\mathit{i},\mathit{o}\mathit{t}\mathit{h}\mathit{e}\mathit{r}\mathit{w}\mathit{i}\mathit{s}\mathit{e}}\end{array}\right.$$

(18)

Here $\mathcal{N}\left(\mathrm{0,1}\right)$is a standard normal random variable generated independently for each coordinate, ⊙ denotes the coordinate-wise (elementwise) product, $\left|x^{*}-x_i\right|$is the per-coordinate distance from $x_i$to the best solution found so far $x^{*}$, so that individuals near $x^{*}$take small steps and distant individuals take larger ones, and ${\eta }_A= 0.25 w + 0.05$ sets the overall magnitude. The candidate $z_i$is accepted only if it improves the objective.

Scaling the step by the distance to the best gives the operator a useful property: individuals already close to $x^{*}$ take small steps and fine-tune, while distant individuals take larger ones and do not get stuck. The step size needs no separate schedule because the population supplies it.

3.9. Maghrib — Greedy Attraction

Maghrib is the shortest window, a boundary event at sunset. It is the most directly exploitative operator: a random 30% of the population move a fraction of the way toward the best solution found so far, keeping the move only if it helps:

$${\mathit{z}}_{\mathit{i}}={\mathit{x}}_{\mathit{i}}+{\mathit{\eta }}_{\mathit{M}}\times \mathit{r}\times \left({\mathit{x}}^{\mathit{*}}-{\mathit{x}}_{\mathit{i}}\right),{\mathit{x}}_{\mathit{i}}=\left\{\begin{array}{c}{\mathit{z}}_{\mathit{i}}, if f\left({\mathit{z}}_{\mathit{i}}\right)<f\left({\mathit{x}}_{\mathit{i}}\right)\\ {\mathit{x}}_{\mathit{i}}, \mathit{o}\mathit{t}\mathit{h}\mathit{e}\mathit{r}\mathit{w}\mathit{i}\mathit{s}\mathit{e}\end{array}\right.$$

(19)

Here $x^{*}$is the best solution to date, $r$ is a uniform random number in [0, 1], and ${\eta }_M = 0.35 (1 - w) + 0.05$ increases as $w$ decreases, so that the attraction is strongest near the end of the run. The candidate $z_i$is accepted only if it improves the objective.

Because ${\eta }_M$is built from $1 - w$rather than $w$, this is the one operator that grows stronger as the run proceeds, which is the conventional behaviour and is placed here deliberately: it is the last intensification before the exploratory Isha phase takes over.

3.10. Isha — Lévy Restart

Isha is the longest window, covering about the final third of the run, and its role is recovery of diversity. Any individual whose stagnation counter has reached the threshold $smax$ = $max(15, 0.1 {Max}_{Iter})$ is displaced by a heavy-tailed Lévy step scaled to the width of the box:

$${\mathit{z}}_{\mathit{i}}={\mathit{x}}_{\mathit{i}}+\mathbf{0.05}\times \mathit{L}\left(\mathit{\beta }\right)\odot \left(\mathit{u}\mathit{b}-\mathit{l}\mathit{b}\right),\mathit{\beta }=\mathbf{1.5}$$

(20)

Here $L\left(\beta \right)$ is a Lévy-flight step—predominantly small but occasionally very large—with stability index $\beta = 1.5$, generated by Mantegna’s method; the coordinate-wise product with ($ub - lb$) scales the step to the width of the search box; and the factor 0.05 keeps its average magnitude moderate.

An individual that has stopped improving therefore usually stays near where it was and sometimes leaves completely. If no individual has stopped improving, the worst 5% are handled instead, so the operator always does something. The move is kept only on improvement.

Placing an exploratory operator last is the point at which PAO departs most clearly from a usual shrinking-only schedule. A monotonically decreasing coefficient cannot re-expand the search once diversity has been lost; a schedule read off the solar cycle can, because the Sun sets as reliably as it rises.

3.11. Boundary Handling

Every candidate produced by any operator is brought back inside the allowed range by clamping (see section3.4): a coordinate that has left the interval is set to the nearest limit it went past. Clamping is fast and never fails, but it piles points up right on the edge, which is itself a known source of built-in search bias. Reflection back into the interior is the usual alternative and is worth testing as part of the sensitivity study.

3.12. Exploration, Exploitation, and the Balance Between them

Table 9 classifies the operators. Two features of the ordering matter. First, the phases run exploration, neutral, balancing, exploitation, exploitation, exploration across the cycle, so the sequence is not monotone. Second, the exploration phase is the last and the longest, which means diversity lost during Asr and Maghrib can be regained before the run ends.

To report the balance with numbers rather than by simply claiming it, population diversity should be recorded at every iteration and turned into exploration and exploitation percentages using the standard measure: the diversity is the average over the dimensions of the mean distance of the individuals from the population median, the exploration percentage is that diversity as a fraction of its largest value seen so far, and the exploitation percentage is the remainder.

Plotting these two curves against the iteration index, with the five windows shaded as in Figure 2, gives direct evidence for or against the claim that the schedule is non-monotone. This is the figure to include in Section 5.7, and it is the one a doubtful reviewer will look for first.

3.13. Control Parameters and Sensitivity

Table 10 lists every quantity that a user could change. They fall into three groups: the two that any population method has, the engine parameters inherited from the base algorithm, and the parameters introduced by PAO itself.

Three sensitivity experiments are needed, and each answers a different question. A one-factor-at-a-time test of values around the defaults, thirty runs per setting, establishes which parameters matter. A test that removes one phase operator at a time (an ablation study), replacing it with the base engine alone, establishes whether each operator is worth keeping; without this the paper cannot show that the five-phase structure contributes anything beyond the base engine underneath. Testing a range of $\phi$ and n values shows whether the schedule family matters at all, and it is the experiment most specific to this algorithm: if performance does not change with φ and n, the astronomy is just for show and would have to be justified some other way.

One point should be given upped rather than left for a referee to raise. Counted fairly, PAO carries more than five primary control parameters, which places it in the highest-parameter category of the parameter-count taxonomy used in Section 1. Fewer parameters are preferable for industrial use, and the natural response is to consolidate: the four subset fractions could share one value, and the three strength coefficients could be expressed as one coefficient times a fixed per-phase ratio. If the sensitivity study shows that doing so costs little, the consolidated version is the one to present.

3.14. Computational Complexity and Evaluation Accounting

Per iteration the engine performs O(N d) arithmetic for mutation and crossover and O(N log N) for the ranking that supplies the $p_{best}$set; each extra operator adds O(k d) for the k individuals it touches. The solar schedule is computed once. Over a whole run,

$$\mathit{O}.\left({\mathit{M}\mathit{a}\mathit{x}}_{\mathit{I}\mathit{t}\mathit{e}\mathit{r}}\times \mathit{N}\times \left(\mathit{d}+\mathbf{log}\mathit{N}\right)\right)$$

(21)

Equation (21) states the time complexity of a single run. Here O$(·)$ denotes asymptotic order, ${Max}_{Iter}$ is the number of iterations, N is the population size, $d$ is the number of variables, and $log N$ accounts for the ranking required to identify the best individuals. The total cost therefore is scaled with the product of the iteration count, the population size, and the sum of the dimension and a small ranking term; it grows linearly in both population size and dimension, as is typical for this class of algorithms.

Objective evaluations are the quantity that must be controlled for a fair comparison, and they are not constant per iteration. The base engine spends $N$, and the active phase spends the amount listed in the last column of Table 9. Scaling those costs by the window shares of Table 4 gives an expected:

$$\mathit{E}.\left[\mathit{F}\mathit{E} \mathit{p}\mathit{e}\mathit{r} \mathit{i}\mathit{t}\mathit{e}\mathit{r}\mathit{a}\mathit{t}\mathit{i}\mathit{o}\mathit{n}\right]\approx 1.12\times \mathit{N}$$

(22)

Equation (22) accounts for the number of objective-function evaluations performed per iteration, which exceeds $N$ because the phase operators introduce additional evaluations. Here $E[·]$ denotes the expected (average) value, $FE$ denotes function evaluations, and N is the population size. The base search performs $N$ evaluations per iteration, and the phase operators add a further amount, yielding approximately 1.12 $N$ on average, or about 12% more. The figure is obtained by weighing each operator’s cost by the frequency with which its phase is active (Table 4); it is the reason a fair comparison must terminate on an evaluation budget rather than on an iteration count.

So, each iteration costs about 12% more than N evaluations. To keep the comparison fair, the corrected code stops on a fixed evaluation budget $\mathit{M}\mathit{a}\mathit{x}\mathit{F}\mathit{E}\mathit{s} = \mathit{N}$ · ${\mathit{M}\mathit{a}\mathit{x}}_{\mathit{I}\mathit{t}\mathit{e}\mathit{r}}$ rather than on the iteration count, so PAO spends the same number of evaluations as its competitors and simply runs proportionally fewer iterations.

3.15. Pseudocode

Figure 5. illustrate the Pseudocode of Pray Algorithm Optimization (PAO).

Every step map onto one box of the flowchart in Figure 3. The switch at lines 21 to 27 is the six-way branch; note that exactly one branch runs per iteration, so the operators never compound within a single generation.

Figure 8. Flowchart of PAO. The solar schedule is computed once before the loop because it depends only on φ and n. Inside the loop the base engine runs every time, without any test and the phase selector activates at most one of the five operators, or none during the forenoon window.

Figure 8. Flowchart of PAO. The solar schedule is computed once before the loop because it depends only on φ and n. Inside the loop the base engine runs every time, without any test and the phase selector activates at most one of the five operators, or none during the forenoon window.

3.16. Algorithmic Novelty Related to PSO, DE, GWO, SCA

Because PAO’s lower layer is an established adaptive base engine (Section 3.2), the case for novelty rests on the astronomical controller and the five phase operators it activates, not on the search engine underneath. This subsection makes that case explicit by contrasting PAO with four widely used metaheuristics — particle swarm optimization (PSO), differential evolution (DE), the grey wolf optimizer (GWO), and the sine–cosine algorithm (SCA) — each chosen because it embodies one of the mechanisms PAO might be mistaken for. The comparison is organized around the single feature that separates PAO from all four: how the balance between exploration and exploitation is scheduled over a run.

PSO moves every particle by one velocity rule that blends inertia with attraction toward the particle’s own best position and the swarm’s global best (Kennedy & Eberhart, 1995). The rule is homogeneous — the same update acts on every particle at every iteration — and its exploration–exploitation balance is governed by an inertia weight that is almost always driven down a fixed, monotone schedule. PAO differs on both counts. It does not apply one rule uniformly: it selects, from five operators, the one suited to the current solar window, and it leaves a forenoon interval in which no extra operator runs at all. And its balance is not a monotone decay but a non-monotone solar profile in which an exploratory phase (Isha) returns at the end of the run. The one point of contact — Maghrib’s greedy attraction toward the incumbent best — resembles PSO’s social pull, but in PAO it acts on a random subset during a single short phase rather than on the whole population throughout.

DE searches by mutating each solution with a scaled difference between population members and then keeping the better of the parent and the trial; its exploration–exploitation balance is carried by a scale factor and a crossover rate, which in adaptive versions are tuned from recent successes, while any population-size reduction it applies decreases monotonically (Storn & Price, 1997). The contrast with PAO is the same in kind as for the other three comparators. DE applies one mutation-and-crossover rule uniformly across the whole run, whereas PAO overlays an external controller that (i) rescales the engine’s step size and mixing rate by phase-specific factors drawn from the solar schedule rather than from search history alone, and (ii) activates, in different windows, five operators that a bare mutation-and-crossover scheme does not contain: opposition-based re-seeding, a pull toward the population median, distance-scaled Gaussian refinement, greedy attraction to the best, and a heavy-tailed Lévy restart of stagnated individuals. Because that operator layer is what is switched on and off, an ablation that disables it returns PAO to its base engine, which is how the paper isolates the contribution.

GWO steers each search agent toward the three highest-ranked agents, with a coefficient vector that decreases linearly from two to zero across the run; this single decay is the whole of its exploration-to-exploitation mechanism (Mirjalili et al., 2014). The contrast is the same in kind as with PSO but sharper: GWO’s control is one hand-set, monotone linear schedule, whereas PAO’s is a non-monotone schedule computed from solar geometry with no fitted constant. GWO applies leadership-following on every iteration; PAO applies qualitatively different operators in different windows. The resemblance is again narrow — Maghrib’s attraction to the single incumbent best echoes GWO’s pull toward its leaders — but it is confined to one phase rather than being the algorithm’s sole engine.

SCA is the most instructive comparison, because it too drives its search with a smooth periodic function: positions are updated through sine and cosine terms whose amplitude, set by a parameter r₁, decreases linearly toward zero over the run (Mirjalili, 2016). Superficially this is very close to PAO’s use of a trigonometric solar signal. Three differences separate them. First, SCA’s amplitude decays monotonically, so it’s swinging up and down shrinks and the method becomes purely exploitative late in the run, whereas PAO’s driving quantity — the shadow ratio, which behaves like cot α — is non-monotone and grows again toward the ends of the cycle, so exploration recovers. Second, SCA’s periodicity is an abstract oscillation whose range depends on a free coefficient that must be chosen, while PAO’s periodicity is the real solar cycle, fixed completely by latitude and day of year, so the schedule carries no constant for the designer to tune. Third, SCA applies one homogeneous sine–cosine update throughout, whereas PAO uses its periodic signal only to partition the run into phases, each of which then runs its own operator.

Two clarifications follow from these comparisons, and both matter how the contribution should be judged. The individual operators PAO switches among are not themselves new — opposition-based learning, Lévy-flight restarts, Gaussian local search, and greedy selection all predate this work — and neither is the base engine. The novelty lies in their organization: a real astronomical signal, parameterized by nothing but a latitude and a date, is used as an external control schedule that is (a) non-monotone, so that exploration can be regained late rather than only lost; (b) phase-structured, so that distinct operators act in distinct, geometrically defined windows rather than a single rule acting uniformly; and (c) free of any fitted schedule constant, so that changing φ or n produces a whole family of valid schedules at no tuning cost. None of PSO, DE, GWO, or SCA has any of these three properties; their control profiles are all single, monotone, hand-set decays. This is also the sense in which PAO is intended to be strongly rather than weakly inspired, in the terms of Section 1.4: the solar geometry is computed from its governing equations and actively sets the search behavior, not merely named after the fact (Aranha et al., 2022; Sörensen, 2015). Table 11 summarizes the comparison of the four axes that define PAO’s contribution.

4. Experimental Setup

This section describes how PAO was tested, so that the study in Section 5 can be reproduced in full. It specifies the four collections of benchmark problems (see section 4.1), the 42 competing algorithms and the control-parameter values taken from their original publications (see section 4.2), and the protocol that keeps the comparison fair — an identical population size, an identical budget of function evaluations, and 30 independent runs for every algorithm on every problem (see section 4.3). It then defines the measurements used to summarize the raw results and the rule for ranking tied values (see section 4.4), the statistical tests used to decide whether a difference between two algorithms is real rather than chance (see section 4.5), and the computing environment (see section 4.6). The design follows the reproducibility conventions now standard in the metaheuristic literature (Derrac et al., 2011; LaTorre et al., 2021), so that every claim made later can be re-derived from the settings given here.

4.1. Benchmark Suites

PAO is evaluated on four complementary collections — one classical and three from the IEEE CEC competition series — chosen so that the test set is neither dominated by easy separable problems nor by any single generation of difficulty.

Classical set (F1–F23): The 23 functions of (Yao et al., 1999), used throughout the swarm-intelligence literature, are included for comparability with earlier work. They fall into three groups: F1–F7 are unimodal (a single optimum) and probe exploitation and convergence speed; F8–F13 are high-dimensional multimodal (a number of local minima that grows with dimension) and probe exploration and the ability to escape local optima; F14–F23 are fixed-, low-dimensional multimodal (2–6 variables, few minima) and probe the exploration–exploitation balance on small landscapes. F1–F13 are scalable and are run at dimension D = 30 (the conventional setting); F14–F23 keep their native dimensions. The search bounds and the known optimum of each function are listed(Yao et al., 1999) in Table A7.

CEC2017 (Pham et al., 2023): Thirty shifted-and-rotated functions are defined; following the common convention the second function (F2) is omitted because its values are numerically unstable across implementations, leaving 29 functions. They comprise 2 unimodal (F1, F3), 7 simple multimodal (F4–F10), 10 hybrid (F11–F20) and 10 composition (F21–F30) problems on the box [−100, 100] in each of D dimensions, with known optima ${f*}_i= 100 i$. The suite is run at D ∈ {10, 30, 50, 100} to expose scaling behavior. The hybrid and composition functions, which glue several rotated sub-functions together, are the ones on which purely exploitative methods stop making progress — exactly the regime the non-monotone schedule of Section 3 is meant to address (see Table B7).

CEC2019 — the “100-Digit Challenge” (Hamid et al., 2026; Price et al., 2018): Ten functions (CEC01–CEC10). The first three carry their own dimensions (CEC01: 9-D, CEC02: 16-D, CEC03: 18-D); CEC04–CEC10 are all 10-D. Every function is shifted and rotated and has a global minimum value of 1, so results are reported as the gap above 1. The suite is deliberately small in dimension but severely ill-conditioned, which stresses fine-grained exploitation late in a run (see Table C7).

CEC2022 ((Salgotra et al., 2025)): Twelve functions on [−100, 100] in each of D dimensions at D ∈ {10, 20}: 1 unimodal (F1), 4 basic multimodal (F2–F5), 3 hybrid (F6–F8) and 4 compositions (F9–F12). As the most recent suite, it provides an up-to-date and widely reported benchmark (see Table D7).

Together the four suites span separable and non-separable, unimodal and multimodal, low- and high-dimensional, and shifted/rotated/hybrid/composition landscapes, giving a total of well over seventy distinct problem instances once dimensions are counted separately.

4.2. Comparator Algorithms, Parameters, Complexity and Budget

PAO is compared against the 42 algorithms of the taxonomy in Figure 1 (listed in Table 2), which together represent the evolutionary, swarm-intelligence, physical-law-based and miscellaneous families. Every comparator is run with the control-parameter values recommended in its original publication; none is re-tuned for the present problems. Re-tuning only PAO, or only its competitors, is the most common source of biased comparisons, so all algorithms are used “as published” and PAO itself is run with the single default configuration of Table 10 on every problem, with no per-suite tuning.

4.3. Fairness Protocol

Three settings are held identical for every algorithm on every problem, so that any observed difference is attributable to the search behavior and not to the resources granted:

(1)

Identical population size. All algorithms use the same N (reported with the results; N = 60 unless a suite’s official protocol prescribes otherwise).

(2)

Identical evaluation budget. Every run terminates when it has consumed MaxFEs = N · ${Max}_{Iter}$ objective-function evaluations, not when it reaches a fixed iteration count. Evaluations are counted individually inside every algorithm, including PAO’s phase operators, which spend extra evaluations per iteration (Section 3.14). Stopping on MaxFEs therefore removes the roughly 12% budget advantage that an iteration-based stop would otherwise give PAO and gives every comparator the same number of chances to improve. The value of MaxFEs used for each suite is reported with the results.

(3)

Thirty independent runs. Each (algorithm, function, dimension) triple is solved in 30 independent runs with independent random seeds, so that the reported statistics reflect the distribution of outcomes rather than a single lucky or unlucky trajectory.

In addition, all algorithms share the same box constraints and the same uniform-random initialization within the feasible region, no algorithm is given the location or value of the optimum, and all bound-handling is by the saturation rule of Section 3.11. The complete list of per-suite N and MaxFEs values is given in Appendix A.

4.4. Performance Metrics and the Tied-Rank Method

For each (algorithm, function, dimension) the 30 final best-so-far objective values are reduced to five statistics: the mean, the standard deviation, the best, the worst, and the median. The mean ± standard deviation is the primary quality/robustness pair; the median, best and worst guard against a mean distorted by a single outlier run and expose skew.

To summarize performance across functions, algorithms are converted to ranks. On each function the $k$ algorithms are ordered from best ($rank 1$) to worst ($rank k$) by their mean value (with the median used to break exact mean ties). Tied results are given the average of the ranks they span — for example, two algorithms tied for positions 3 and 4 both receive rank 3.5 — so that the ranks on every function still sum to $k(k + 1)/2$ and no algorithm is arbitrarily favoured by the ordering of ties. Each algorithm’s average rank over all functions is then the single number used to compare methods and to build the diagrams of section 4.5. This tied-rank convention is the one required by the Friedman procedure and is what makes the rank-based tests below valid.

4.5. Statistical Tests

Raw means can differ by amounts too small to be meaningful, so all conclusions in Section 5 are supported by non-parametric statistics, which make no normality assumption about the 30-run samples (Derrac et al., 2011).

Wilcoxon rank-sum test (pairwise): For every function, PAO’s 30 results are compared with each competitor’s 30 results using the two-sided Wilcoxon rank-sum (Mann–Whitney U) test at the $\alpha = 0.05$ significance level. Each comparison is scored as a win (+), tie (=) or loss (−) for PAO, and the overall +/=/− tally against each algorithm is reported.

Friedman test (multiple comparison): The average ranks of section 4.4 are tested with the Friedman test, whose null hypothesis is that all algorithms have equal average rank. The Friedman statistic, its p-value and the resulting ranking are reported; a small p-value licenses the post-hoc analysis.

Post-hoc analysis: When Friedman rejects the null, a Holm step-down procedure (with PAO as the control method) identifies which competitors PAO significantly out- or under-performs while controlling the family-wise error rate across the many pairwise comparisons; the adjusted p-values are reported.

Critical-difference (CD) diagrams: Following (Demšar, 2006) average ranks are drawn on a single axis with a critical difference bar, )CD =$q\alpha · \surd ( k(k + 1) / (6 N꜀\left)\right),$ where k is the number of algorithms and N꜀ the number of functions. Algorithms whose average ranks differ by less than CD are joined by a bar and are statistically indistinguishable; those not joined differ significantly. The CD diagram gives an at-a-glance summary of the whole comparison and is the figure a reviewer is most likely to inspect first .

Where a suite is run at several dimensions, the tests are applied separately per dimension and then on the pooled ranks, so that both dimension-specific and overall conclusions are supported.

4.6. Computing Environment

All algorithms, including PAO, are implemented in MATLAB [R2020a] and run on the same machine — processor: Intel Core i7-10700 @ 2.90 GHz, 8 cores; Memory: 32 GB RAM; operating system: Microsoft Windows 11 pro — under identical conditions, with each run executed on a single core and no parallelism across runs, so that the wall-clock times used to validate the complexity analysis of section 3.14 are comparable. Reported run-times are the average over the 30 runs. The PAO source code and the scripts that regenerate every table and figure are available at [ https://www.mathworks.com/matlabcentral/fileexchange/184442-pao-for-global-and-constrained-engineering-optimization ].

5. Results and Discussion on Benchmark Functions

This section reports how PAO performed on the four collections of standard test problems. For each collection it gives the average result of every algorithm, ranks them, and reports whether PAO’s advantage is statistically significant. It then looks at the results from several angles: convergence curves, boxplots of run-to-run consistency, the shifting balance between searching and refining, and tests of how the method behaves when the problem gets larger, when the population size changes, and how much time it takes. Where PAO loses, the section says so and explains why.

5.1. Classical Benchmark Functions_23 Functions

This section evaluates the Pray Algorithm Optimization (PAO) against a large contemporary pool of forty-three metaheuristics on the classical set of twenty-three benchmark functions (F1–F23). Performance on every function is the mean fitness over the 30 runs. Per-function ranks were then assigned in ascending order of mean fitness (rank 1 = best), with ties resolved by the standard fractional (mid-rank) convention, and each algorithm’s overall standing was taken as its mean rank across the 23 functions. All supporting numbers are collected in Appendix A.

Averaged over the 23 functions, PAO attains a mean rank of 13.152, placing it second of the 43 algorithms. It is bettered only by SOA, whose mean rank of 11.261 heads the field, and it leads a close chasing group formed by WAA (13.870, 3rd) and DHOLE (14.630, 4th). The Friedman overall test confirms that these rank differences are not attributable to chance (χ² = 311.91, $df = 42, p < 0.001)$, so the ordering in Table A1 is statistically meaningful rather than an artefact of sampling. PAO can therefore be described, without embellishment, as a strong second-best optimizer on this suite — highly competitive with the leader but not its equal on the combined measure.

The Wilcoxon rank-sum test (α = 0.05) applied function-by-function paints a consistent picture. Added together over all 42 competitors and 23 functions, PAO records 641 wins, 42 ties and 283 losses — a 66.4%-win rate (Table A4). Resolved competitor by competitor, PAO holds a winning record against 33 of the 42 competitors, an even record against 2 (AO and EMRA), and a losing record against only 7. Those seven are SOA and CSA (both 9–14 against PAO), followed by GOA, DHOLE, TSO, SO and HBA, each of which edges PAO by only one or two functions. In other words, exactly one algorithm (SOA) both out-ranks PAO on the mean-rank metric and beats it head-to-head by a clear margin; PAO’s runner-up position is thus robust to the choice of statistical lens.

PAO’s strongest placements are the harder fixed-dimension multimodal landscapes, where it ranks first or near-first on F18, F21, F22 and F23; the per-function ranks are given in Appendix A.

On the easy functions a large fraction of the 43 algorithms drive the mean fitness exactly onto the global optimum — that is, all 30 of their runs converge to the optimal value — so their mean finesses are numerically identical. Under the fractional-ranking rule, a block of k algorithms sharing the best mean does not occupy ranks 1 through k; every member is instead assigned the single averaged rank (k+1)/2. The optimum therefore collapses a whole span of ranks into one value. The effect is severe on this suite (Table A2): the top block contains 18 algorithms on F11 (each assigned rank 9.5), 17 on F9 and again 17 on F16 (rank 9.0), 16 on F10 (rank 8.5), 15 on F19 (rank 8.0), 14 on F17 (rank 7.5), 10 on F14 (rank 5.5) and 8 on each of F1–F4 (rank 4.5). A direct consequence is a loss of resolution in the mean-rank metric itself: on F11, for example, the ranking cannot distinguish the best optimizer from the eighteenth-best, because all eighteen are recorded at 9.5. The between algorithm variance on such functions is largely erased, and the aggregate ranking is decided disproportionately by the minority of functions on which ties do not form.

The compression acts on PAO in both directions: on several separable functions its mean sits a hair above the tie block and is displaced below algorithms it effectively matches, while elsewhere the same rule works in its favor. Appendix A gives the per-function figures.

Because so many ranks depend on the optimum, PAO’s aggregate score is convention dependent. Recomputing the mean rank under the optimistic (minimum) and pessimistic (maximum) tie rules bounds PAO between 11.913 and 14.391 (Table A1). Under the pessimistic rule the reported 13.152 is preserved and PAO remains second; under the optimistic rule PAO would be narrowly overtaken by DHOLE (whose min-tie score of 11.826 benefits from its membership in many optimum blocks), dropping PAO to third. SOA, by contrast, remains first under every convention (min-tie 10.130). The second-place finish is thus stable under the standard fractional rule used throughout this study, while the exact second/third boundary is acknowledged to be tie-rule sensitive.

Taken together, PAO’s mean rank of 13.152 is best read as a conservative estimate of its quality: it is lowered by tie compression precisely on the functions PAO solves to optimality and lifted toward first only by its tie-free, outright wins on the harder multimodal problems. SOA edges ahead because it exploits both fronts at once — it wins the separable unimodal functions clearly (its F1 mean reaches 3.89 × 10−²²³) and simultaneously sits inside the multimodal optimum blocks — whereas PAO trades block membership on the easy functions for dominance on the difficult ones. The same tie structure explains a mild non-transitivity in the results: PAO out-ranks GOA, TSO and HBA on the mean-rank metric yet is slightly net-negative against them head-to-head, because averaged ranks and pairwise counts weight the compressed functions differently. The conclusion holds: on the classical 23-function suite PAO is a robust second, its comparative advantage concentrated on multimodal and composite optimization.

5.2. CEC2017 Benchmark Functions_30 Functions

The CEC-2017 suite was run at dimension D = 30 for the same pool of forty-three algorithms, each for 30 independent runs of 500 iterations. Following the standard convention adopted in Section 4.1, the second function (F2) is excluded because its shifted-and-rotated form is numerically unstable across implementations, leaving twenty-nine functions: two unimodal (F1, F3), seven simple multimodal (F4–F10), ten hybrid (F11–F20) and ten compositions (F21–F30). Performance on each function is the solution error f(x) − f* over the 30 runs (the known optima are f*ᵢ = 100i, so lower error is better and zero denotes the optimum located). Ranking follows the convention of Section 5.1; supporting numbers are in Appendix B.

PAO ranks first. Averaged over the 29 functions, PAO attains a mean rank of 2.172, placing it first of the 43 algorithms clearly. The nearest competitors trail well behind — DOA (5.655, 2nd), NDO (5.966, 3rd), PGA (7.034, 4th), EHO (8.586, 5th) and GOA (9.207, 6th) — so PAO’s average rank is better than the runner-up’s by more than three full rank positions. The dominance is repeated in outright wins: PAO ranks first on 19 of the 29 functions, whereas no other algorithm is first on more than three. The Friedman overall test rejects rank equality overwhelmingly (χ² = 949.83, df = 42, p ≈ 8 × 10^⁻¹⁷²), so this ordering is not a chance result of sampling. On this suite PAO can be described, without qualification, as the strongest optimizer in the pool.

The two-sided Wilcoxon rank-sum test (α = 0.05), applied function-by-function, tells the same story. Aggregated over all 42 opponents and 29 functions, PAO records 1150 wins, 56 ties and 12 losses — a 94.4%-win rate (Table B4). Resolved opponent-by-opponent, PAO holds a net-positive record against every one of the 42 competitors without exception and sweeps 16 of them 29–0–0. Even the strongest rivals are beaten decisively: DOA, NDO, EHO, PGA and PSO — the only algorithms to take a single function from PAO — still concede it a net advantage of between +20 and +23. In short, no algorithm in the pool out-ranks PAO or holds an even head-to-head record against it; the first-place finish is robust to the choice of statistical lens.

The advantage is broad rather than concentrated: PAO ranks first on both unimodal functions and on six of the seven simple-multimodal functions and leads the hybrid and composition groups. Per-function errors are listed in Appendix B.

Its one clear weakness is F10 (Shifted-and-Rotated Schwefel, rank 14), whose deceptive geometry places the global basin far from the strong local basins and traps PAO’s late refinement; see Appendix B.

The first-place finish here stands in deliberate contrast to PAO’s second place on the Classical-23 suite (Section 5.1), and the two results are consistent once the difference between the suites is understood. On the easy, separable Classical functions many algorithms drive the mean exactly onto the optimum, so their ranks collapse into large ties and the metric loses resolution — an effect that masked PAO’s true accuracy there. The CEC-2017 functions are shifted, rotated, hybridized and composed specifically to defeat this: almost no algorithm reaches the exact optimum, ties at the optimum vanish, and the ranking therefore has full resolution. Under that sharper lens PAO’s genuine search ability is expressed directly. Tellingly, the algorithms that led the Classical suite by sitting inside optimum tie-blocks collapse here — SOA falls from 1st on Classical-23 to 32nd, and CSA to 42nd of 43 — which indicates their Classical lead reflected tie-compression rather than search power. On the harder, more discriminating benchmark, PAO is the clear leader.

Taken together, the mean rank of 2.172, the 19 first-place finishes, the 94.4% pairwise win rate and the net-positive record against all 42 competitors establish PAO as the best-performing optimizer on CEC-2017 at D = 30. Its comparative advantage is concentrated exactly where it was designed to be strong — the multimodal, hybrid and composition landscapes — and its sole material weakness (the Schwefel function F10) is specific and reproducible. Runtime cost is discussed in Section 5.11; on the quality metrics that matter, PAO leads the field.

5.3. CEC2019 Benchmark Functions_10 Functions

The CEC-2019 “100-Digit Challenge” comprises ten functions (CEC01–CEC10) evaluated for the same pool of forty-three algorithms, each for 30 independent runs of 500 iterations. The first three are real-world problems carrying their own dimensions — Storn’s Chebyshev polynomial fitting (9-D), the inverse Hilbert matrix (16-D) and the Lennard-Jones minimum-energy cluster (18-D) — while CEC04–CEC10 are shifted-and-rotated 10-D problems (Rastrigin, Griewank, Weierstrass, modified Schwefel, expanded Schaffer F6, Happy Cat and Ackley). Every function has a global minimum value of 1, so performance is reported as error (the objective gap above 1), for which lower is better. The suite is small in dimension but severely ill-conditioned, which stresses fine-grained exploitation late in a run. Ranking follows Section 5.1; supporting numbers are in Appendix C.

PAO ranks first. Averaged over the ten functions, PAO attains a mean rank of 4.55, placing it first of the 43 algorithms — and by a striking margin, since the runner-up NDO sits at 10.05, more than double PAO’s figure. The next places are tightly grouped closely and far behind: PGA (10.20, 3rd), DOA (10.80, 4th), PSO (10.85, 5th) and GOA (11.50, 6th). The Friedman overall test rejects equal ranking decisively (χ² = 228.86, df = 42, p ≈ 1.5 × 10^⁻²⁷). PAO wins here through consistency, not outright wins: first on only three functions (F4, F5, F9) but top-six on eight of ten, while no competitor is uniformly strong — the chasing group clusters near an average rank of ten precisely because each of them is weak on some subset of the functions. PAO’s low average rank is therefore earned by rarely placing poorly.

The two-sided Wilcoxon rank-sum test (α = 0.05), applied function-by-function, confirms the picture. Aggregated over all 42 opponents and ten functions, PAO records 362 wins, 31 ties and 27 losses — an 86.2%-win rate (Table C4) — and holds a net-positive record against every one of the 42 competitors without exception, sweeping 18 of them 10–0–0. The closest rivals still concede ground: HBA, PGA and NDO, the strongest of the field, each finish only +4 in net terms against PAO. No algorithm out-ranks PAO or holds an even head-to-head record against it, so the first-place finish is robust to the choice of statistical lens.

The advantage is concentrated on the shifted-and-rotated 10-D functions: PAO ranks first on Rastrigin, Griewank and Happy Cat, beating all 42 competitors on each. Appendix C gives the per-function detail.

The exceptions are more numerous here than on CEC-2017. PAO is weakest on F6 (Weierstrass, rank 11), whose fractal landscape defeats late-stage refinement, and only mid-pack on the three real-world functions. Appendix C lists them.

The first-place finish reinforces the CEC-2017 result (Section 5.2): on hard, well-conditioned benchmarks that discriminate sharply between algorithms, PAO leads the field. The difference in character — a consistency-driven win here versus a win-almost-everything result on CEC-2017 — reflects the nature of the 100-Digit suite, where the severe conditioning of individual functions means no single algorithm dominates all ten, and the best overall method is the one that is strong across the board. By that measure PAO is the clear winner, just as it is on CEC-2017; the contrast with its masked second place on the easy Classical-23 suite (Section 5.1) is again explained by the absence here of the optimum-tie compression that flattened the Classical ranking.

Taken together, the mean rank of 4.55 — less than half the runner-up’s — the 86.2% pairwise win rate and the net-positive record against all 42 competitors establish PAO as the best-performing optimizer on CEC-2019. Its strength lies on the shifted-and-rotated multimodal functions, and its limitations are specific: the Weierstrass function and the specialized real-world problems. Runtime is discussed in Section 5.11; on the quality metrics that matter, PAO again leads.

5.4. CEC2022 Benchmark Functions_12 Functions

The CEC-2022 suite — the most recent of the four — comprises twelve functions on [−100, 100]: one unimodal (F1, shifted-rotated Zakharov), four basic multimodal (F2–F5: Rosenbrock, expanded Schaffer F6, non-continuous Rastrigin, Levy), three hybrid (F6–F8) and four composition (F9–F12), with known optima f*ᵢ ranging from 300 to 2700. The suite is defined at D ∈ {10, 20}; the results reported here are for D = 10. Performance on each function is the mean solution error f(x) − f* over the runs (lower is better; zero denotes the optimum located), and ranking follows Section 5.1. Supporting numbers are in Appendix D.

PAO ranks first, but narrowly. Averaged over the twelve functions, PAO attains a mean rank of 4.375, the best in the pool — but here the margin is slim. GOA follows at 4.667 and APO at 5.833, so PAO, GOA and APO form a clearly separated leading trio ahead of PGA (8.708) and DOA (9.083). The Friedman overall test confirms that the overall ranking is highly significant (χ² = 378.08, df = 42, p ≈ 1.2 × 10^⁻⁵⁵), but the gap between PAO and GOA at the very top is small enough that the two are best regarded as sharing the first tier rather than cleanly separated. PAO takes outright first place on four functions (F1, F4, F5, F6); GOA, APO and DOA each win two.

The two-sided Wilcoxon rank-sum test (α = 0.05) reinforces this reading. Aggregated over all 42 opponents and twelve functions, PAO records 445 wins, 32 ties and 27 losses — an 88.3% win rate (Table D4) — and holds a net-positive record against 41 of the 42 competitors, sweeping 20 of them 12–0–0. The single exception is instructive: against GOA, PAO is exactly level — five wins, two ties and five losses — so although PAO edges GOA on the average rank metric, the two are statistically indistinguishable head-to-head. APO, the third member of the leading trio, is beaten by a clear +4 (7–2–3). PAO is therefore the top-ranked method on this suite, but honestly it shares the summit with GOA rather than standing alone.

PAO is strongest on the unimodal, basic-multimodal and first-hybrid functions, solving Zakharov and Levy to the optimum on every run. Appendix D gives the per-function results.

The exceptions cluster on two families, the worst being Rosenbrock (rank 11), where a narrow curved valley leaves a mean error of 6.81 although the best run reaches the optimum — a consistency shortfall rather than an inability to solve it. Appendix D has the detail.

This is PAO’s third consecutive first-place finish on a modern, discriminating benchmark, after CEC-2017 and CEC-2019, and it reinforces the overall conclusion that PAO is a leading optimizer on hard problems. The narrower margin here — a near-tie with GOA rather than the wide lead seen on CEC-2017 — reflects both the strength of GOA on this particular suite and that CEC-2022 is the newest and, at the top of the field, the most closely contested collection. As on the other CEC suites, the absence of optimum-tie compression means the ranking has full resolution, so PAO’s first place is a genuine reflection of accuracy rather than the artefact that flattened its Classical-23 standing (Section 5.1).

Taken together, the best mean rank in the pool (4.375), the 88.3% pairwise win rate and the net-positive record against 41 of 42 competitors place PAO at the top of the CEC-2022 ranking, with GOA as a statistically level co-leader. PAO’s advantage lies on the unimodal, basic multimodal and hybrid functions — several solved to the optimum on every run — while its limitations are confined to Rosenbrock and two composition functions. Runtime is discussed in Section 5.11; on quality, PAO leads this suite too, if narrowly.

5.5. Convergence and Distribution of Results

5.5.1. Classical-23

Figure 6 traces the mean best-so-far fitness of every algorithm on each of the 23 functions, with PAO highlighted against the field and a logarithmic vertical axis wherever fitness is strictly positive. The curves make PAO’s search dynamics explicit and complement the end-point ranking of Section 5.1.

On the unimodal and scalable functions PAO exhibits a clean, multi-phase descent. An exploratory phase over the first 10–20% of the budget cuts the error by orders of magnitude — on F1 the mean best-so-far falls from 5.10 × 10⁴ at iteration 1 to about 1.75 by iteration 100 — followed by a sustained, near log-linear refinement phase and a late plateau at the optimum (2.97 × 10−¹² by iteration 500). The same signature appears on the steep F5 valley, where the curve descends from 1.80 × 10⁸ to 32.4, and on the highly multimodal F9, where it falls from 4.03 × 10² to 0.895. No oscillation, stagnation spike or divergence is observed on any PAO trace.

On the fixed-dimension multimodal functions PAO locates and settles into the global basin within a few tens of iterations, consistent with its top ranks there: on F23 the curve moves from −1.28 to the optimum −10.54 by roughly iteration 50, with F18, F21 and F22 behaving similarly. The convergence-iteration summary (Table A5) quantifies this: PAO stabilizes very early on the fixed-dimension set (iteration 31 on F16, 25 on F18) yet continues to improve late on the hardest separable functions (iteration 490 on F4, 497 on F8). This combination — early lock-on where the landscape permits, patient progress where it does not — indicates genuine iterative optimization rather than premature convergence.

The strongest exploiters (SOA, CSA) sit on the numerical floor for most of the run. On the separable functions their best-so-far trace hits the floor within the first logged iterations — SOA’s trace on F1, F5 and F9 is already pinned at zero from iteration 1, although its terminal mean on F1 is 3.89 × 10−²²³ — so an iteration-by-iteration curve comparison against these algorithms is uninformative on those functions, and their relative standing there is better read from the terminal accuracy in Table A3. On the harder multimodal and fixed-dimension functions (F8, F14–F23), where all traces are well resolved, PAO’s curve is competitive with, and frequently ahead of, the leading algorithms — the visual counterpart of its first-place ranks on those problems.

Figure 6. Convergence behavior on the Classical-23 benchmark suite (30-run mean best-so-far, 500 iterations) PAO highlighted against the field.

Figure 6. Convergence behavior on the Classical-23 benchmark suite (30-run mean best-so-far, 500 iterations) PAO highlighted against the field.

Where the convergence curves describe the typical run, the boxplots in Figure 7 describe the spread across the 30 runs and therefore the robustness of each optimizer. For legibility the panels compare PAO against the eight leading algorithms by mean rank; each box spans the inter-quartile range, the whiskers reach the non-outlier extremes, and individual markers denote outliers, with a logarithmic axis where fitness is positive.

On the functions that drive its ranking — the fixed-dimension multimodal set F16–F23 and the step function F6 — PAO produces tight, low boxes with negligible height, meaning it converges to the same optimum on every run. Its dispersion of these problems matches or improves upon that of the leading algorithms, so PAO’s strong average ranks their rest on repeatable behavior rather than on a few fortunate runs.

On the easy separable functions PAO’s box sits slightly above the leaders and is visibly taller. F9 is the clearest case: PAO’s 30 final values have mean 0.895 and standard deviation 0.990, ranging from a best of 1.85 × 10−⁶ to a worst of 3.98, indicating that a minority of runs settle on a Rastrigin local basin instead of the global optimum. This dispersion is exactly the mechanism identified in Section 5.1: it lifts PAO’s mean just outside the optimum tie block and, through rank compression, drops its rank on those functions. The leaders (SOA, CSA) by contrast show degenerate zero-height boxes on the same problems, every run having reached the optimum.

PAO shows few high outliers, and its whiskers stay close to the box, so even its worst runs are bounded. On most functions the median lies near the center of the box, indicating symmetric run distributions. The boxplots corroborate the ranking of Section 5.1: PAO is highly consistent on multimodal and fixed-dimension problems, and competitively consistent, if slightly more dispersed, on the separable functions where the optimum tie block depresses its rank.

Figure 7. Distribution of final fitness over 30 independent runs Classical-23 suite PAO versus the eight leading optimizers (lower is better).

Figure 7. Distribution of final fitness over 30 independent runs Classical-23 suite PAO versus the eight leading optimizers (lower is better).

5.5.2. CEC-2017

Figure 8 traces the mean best-so-far error of every algorithm on each of the 29 CEC-2017 functions (F2 omitted), with PAO highlighted and a logarithmic vertical axis.

On the unimodal and simple-multimodal functions PAO shows a steep, multi-order-of-magnitude descent. On F1 the mean error falls from 1.05 × 10^¹¹ at iteration 1 to about 1.19 × 10^⁷ by iteration 100 and 1.06 × 10^³ at the end — roughly eight orders of magnitude — and a similar shape appears on F5 (622 → 18.7) and F11 (3.46 × 10^⁴ → 66.9). The trajectory is front-loaded: most of the reduction occurs in the first 100–150 iterations, after which PAO settles into a refinement phase rather than stagnating. No oscillation or divergence is visible on any PAO trace.

On the composition set, PAO descends promptly to the basin floor — F21 moves from 785 to 218 and F23 from 1485 to 368 — and its curve sits at or near the bottom of the field on the large majority of panels, the visual counterpart of its first-place ranks there. The exception is visible too: on F10 the curve descends only modestly (9447 → 4711) and plateaus well above the leaders, which is exactly the rank-14 outcome discussed in Section 5.2. The panels show early, clear progress on almost every function, with the Schwefel function the single case where PAO’s late-stage refinement is out-paced.

Figure 8. Convergence behavior on the CEC-2017 benchmark suite (30-run mean best-so-far error, 500 iterations, D = 30 Colors as in Figure 6. F2 is omitted per the standard convention.

Figure 8. Convergence behavior on the CEC-2017 benchmark suite (30-run mean best-so-far error, 500 iterations, D = 30 Colors as in Figure 6. F2 is omitted per the standard convention.

Where the convergence curves describe the typical run, the boxplots in Figure 9 describe the spread of the final error across the 30 runs. For legibility the panels compare PAO against the seven leading algorithms by mean rank; each box spans the inter-quartile range, the whiskers reach the non-outlier extremes, and a logarithmic axis is used.

PAO is strikingly consistent on the composition set. Its relative dispersion (standard deviation divided by mean error) is only 2.4% on F21, 1.4% on F23, 1.0% on F24 and 0.5% on F25 — tight, low boxes with negligible height, and the lowest in the field on those functions. On the simple-multimodal functions that it solves to near-zero it is equally repeatable, with small spreads on F6 (best 5.30 × 10^⁻³) and F9 (best 1.18 × 10^⁻⁸). PAO’s strong average ranks therefore rest on repeatable behavior rather than on a few fortunate runs.

The dispersion is larger on a few functions. On F1 (Bent Cigar) the box is tall relative to its center — mean error 1063 with a standard deviation of 1187 — although PAO still ranks first there, and on the weak function F10 the box sits high, consistent with its rank-14 outcome. The one composition function with genuinely large spread is F22 (relative dispersion 1.74), where runs occasionally settle in different basins and AOA and APO edge PAO. Elsewhere PAO exhibits few high outliers and its whiskers remain bounded, so its worst runs are rarely catastrophic.

The boxplots corroborate Section 5.2: PAO is not only the most accurate optimizer on this suite on average, but also among the most consistent, with tight low-error distributions across the hybrid and composition families where it leads and only a small number of functions — chiefly F10 and F22 — on which its run-to-run behavior is more variable.

5.5.3. CEC-2019

Figure 10 traces the mean best-so-far error of every algorithm on each of the ten CEC-2019 functions, with PAO highlighted and a logarithmic vertical axis.

On the shifted-and-rotated functions PAO shows a rapid, early, rapid descent. On Rastrigin (F4) the mean error falls from 1.97 × 10^⁴ at iteration 1 to about 43 by iteration 100 and 2.02 at the end, and on Happy Cat (F9) it drops from 4.13 × 10^³ to 2.34, converging to the optimum basin within roughly a hundred iterations; on Griewank (F5) the curve settles to 1.002, the optimum value. Most of the reduction on these functions occurs in the first 100–150 iterations, after which PAO refines rather than stagnates, and no oscillation or divergence is visible on any PAO trace.

Two behaviors visible in the figure match the ranking of Section 5.3. On the real-world problem F1 (Chebyshev) PAO descends about eight orders of magnitude — from 1.96 × 10^¹² to 4.20 × 10^⁴ — yet several competitors reach lower values, so the steep curve still yields only a sixth-place finish. On Weierstrass (F6) the curve descends only modestly (14.3 → 7.83) and plateaus above the leaders, the visual counterpart of its rank-11 outcome. Read together, the panels show decisive early progress on the multimodal functions PAO wins, with the fractal Weierstrass landscape the one case where its late-stage refinement is clearly out-paced.

Where the convergence curves describe the typical run, the boxplots in Figure 11 describe the spread of the final error across the 30 runs. For legibility the panels compare PAO against the seven leading algorithms by mean rank; each box spans the inter-quartile range, the whiskers reach the non-outlier extremes, and a logarithmic (or symmetric log, where values are non-positive) axis is used.

On several functions PAO’s 30 runs are almost identical. On the inverse Hilbert (F2) and Lennard-Jones (F3) problems its standard deviation is below 10^⁻¹⁴ — the box collapses to a line — and on Griewank (F5) and Happy Cat (F9) the relative dispersion is only 0.4% and 0.1% respectively. PAO reliably converges to the same basin on these problems; on F5 and F9 that basin is the best in the field, whereas on F2 and F3 it is a consistently reached but non-winning value, which is why the near-perfect repeatability there still yields a fifth-place rank.

Two functions show real dispersion. On Ackley (F10) the distribution is effectively bimodal: the best run drives the error to 2.34 × 10^⁻⁸ (the optimum) while the worst stalls near 20, so PAO solves the function on some runs and not others — yet its median is still good enough for third place. On the Chebyshev problem (F1) the spread is the widest on the suite, with runs ranging from 3.14 × 10^⁴ to 3.49 × 10^⁵. Weierstrass (F6), by contrast, is tightly distributed (relative dispersion 7%) but sits high, consistent with its weak rank. Rastrigin (F4) shows moderate spread yet remains first.

The boxplots corroborate Section 5.3: PAO is the most accurate optimizer on this suite on average and is highly repeatable on most functions, with only Ackley (bimodal) and the Chebyshev problem (wide) showing substantial run-to-run variation and Weierstrass the one function on which its consistent behavior settles at an uncompetitive value.

5.5.4. CEC-2022

Figure 12 traces the mean best-so-far error of every algorithm on each of the twelve CEC-2022 functions (D = 10), with PAO highlighted and a logarithmic vertical axis.

On the functions it wins PAO shows a fast, deep descent that reaches the optimum. On Zakharov (F1) the mean error falls from 7.77 × 10^⁴ at iteration 1 to the global optimum by the end of the run, and on the first hybrid function (F6) it drops nine orders of magnitude — from 7.63 × 10^⁸ to 0.45 — with most of the reduction complete within the first 100–150 iterations. Levy (F5) is driven to exactly zero on every run, and non-continuous Rastrigin (F4) settles near the optimum (mean error 2.69). No oscillation or divergence is visible on any PAO trace.

The functions on which PAO ranks poorly are equally visible. On Rosenbrock (F2) the curve descends quickly to about 8 and then plateaus — from 1.64 × 10^³ to 6.81 — the flat tail above the optimum being the counterpart of its rank-11 finish. On the composition functions F10 and F11 the curves settle onto raised basins (F10 to an error near 104, F11 to about 210), where several competitors reach lower values. Read together, the panels show PAO reaching the optimum decisively on the unimodal, basic and hybrid functions and stalling above it on Rosenbrock and part of the composition set.

Where the convergence curves describe the typical run, the boxplots in Figure 13 describe the spread of the final error across the 30 runs. For legibility the panels compare PAO against the seven leading algorithms by mean rank; each box spans the inter-quartile range, the whiskers reach the non-outlier extremes, and a logarithmic axis is used.

On the functions PAO wins, its distribution collapses onto the optimum. On Levy (F5) all thirty runs return the exact optimum (zero spread), and on Zakharov (F1) the standard deviation is at the level of numerical noise (10^⁻¹⁴) — the boxes are single lines at the floor, the best in the field. Non-continuous Rastrigin (F4) and the first hybrid (F6) show only modest spread while remaining first. PAO’s winning ranks on these functions therefore rest on almost identical optimal behavior rather than on occasional lucky runs.

patterns explain the weaker ranks. Composition function F11 is effectively bimodal: the best run reaches the optimum (error 0) while others settle near 300, so PAO solves it on some runs and not others. On composition function F10 the picture is different — PAO’s median is strong (fourth in the field) but a few high-outlier runs inflate its mean to a rank-10 value, an asymmetric distribution rather than a uniformly poor one. Rosenbrock (F2) is the one genuinely mid-pack case, with a moderate, consistent spread that sits above the leaders. Elsewhere PAO’s whiskers stay bounded and high outliers are rare.

The boxplots corroborate Section 5.4: PAO is the top-ranked optimizer on this suite and is almost identically optimal on the functions it wins, with its weaker composition results driven by bimodality (F11) or a handful of outlier runs (F10) rather than by broadly poor behavior — and with GOA close enough that the two share the leading position.

5.6. Qualitative Analysis: Search History, Trajectory, Average Fitness, Exploration/Exploitation Percentages

A qualitative look at one representative run explains how PAO searches. Four standard diagnostics are recorded on the classical functions in two-dimensional form, so the population can be seen directly: the search history (the positions visited by all search agents over the whole run), the trajectory (the first coordinate of a single agent across iterations), the average fitness of the whole population across iterations, and the exploration/exploitation percentages derived from population diversity. Figure 11 and Figure 12 show the full five-panel composite for a unimodal and a multimodal function; the remaining functions are collected in Appendix E.

Each composite sets the two-dimensional landscape (panel a) beside the search history (panel b), the trajectory (panel c), the population-average fitness (panel d) and the convergence curve (panel e). The search-history panels show the expected pattern: on the unimodal functions the visited points concentrate tightly around the global optimum, a visible signature of exploitation, whereas on the multimodal functions the agents first scatter across several basins before contracting onto the best one. The trajectory panels show large-amplitude oscillations in the early iterations — the agent sweeping across the domain during exploration — that damp progressively into a small, stable neighborhood as the search settles, the sudden reductions in amplitude reflecting the non-monotone astronomical schedule of Section 3. The average-fitness panels fall and then flatten while tracking the best-so-far curve closely, which confirms that the whole population converges toward the optimum rather than a single lead agent dragging a dispersed swarm. Together the three panels depict a search that explores broadly at first and then commits decisively to refinement (see Figure 14 and Figure 15.).

The balance is quantified from the dimension-wise population diversity: at each iteration the exploration percentage is the current diversity relative to its maximum over the run, and the exploitation percentage is its complement, the two summing to 100%. Averaged over a run, these percentages summarize the search regime. For PAO the balance is strongly exploitation-dominated on most functions — across the five functions studied it averages 15.5% exploration to 84.5% exploitation — which is consistent with the fast early descent seen in the convergence curves of Section 5.5 and with PAO’s strength on functions that reward rapid refinement. The per-function values are given in Table 12 and Figure 16.

The single most informative feature of Table 10 is that PAO’s balance is adaptive, not fixed. On the unimodal Sphere (F1) and the fixed-dimension Kowalik and Hartman functions (F15, F19) it exploits almost throughout (exploration 3–7%), but on the highly multimodal, deceptive Schwefel function (F8) — whose global optimum lies far from the strong local basins — it more than doubles its exploration to 52%, correctly sustaining search where premature convergence would be fatal. Set against the 43-algorithm field the pattern is sharper still: PAO explores far less than the average competitor on the functions that reward exploitation (3.8% vs a field mean of 13.1% on F1, 3.1% vs 32.9% on F19) yet more than the field on F8 (52.3% vs 47.0%). This selective allocation — exploring only where the landscape demands it and refining elsewhere — is the behavioral counterpart of the accuracy reported in Section 5.2, Section 5.3, Section 5.4 and Section 5.5 and follows directly from the astronomically scheduled search. The complete per-function composites and diversity curves are provided in Appendix E.

5.7. Statistical Analysis: Friedman ANOVA, Wilcoxon Summary and Critical-Difference Diagrams

Rankings alone do not establish that one optimizer is better than another; they must be tested. Three complementary tests are applied to the four benchmark suites of Section 5.1, Section 5.2, Section 5.3 and Section 5.4. Friedman ANOVA is an overall test asking whether the 43 algorithms differ at all in their average ranks. The Wilcoxon rank-sum test (α = 0.05) is applied function-by-function to every PAO-versus-competitor pair and added together into wins, ties and losses. Finally, critical-difference (CD) diagrams based on the post-hoc Bonferroni–Dunn procedure, with PAO as the control method, show which algorithms are separated from PAO by a statistically meaningful margin. Table 13 collects the outcomes and Figure 17 presents the CD diagrams; per-suite detail is given in Appendix F.

The overall test is decisive on every suite. The test statistics are χ² = 311.91 on Classical-23 (df = 42, p = 6.35 × 10^⁻⁴³), χ² = 981.42 on CEC-2017 (p = 2.17 × 10^⁻¹⁷⁸), χ² = 228.86 on CEC-2019 (p = 1.48 × 10^⁻²⁷) and χ² = 378.08 on CEC-2022 (p = 1.24 × 10^⁻⁵⁵). Every p-value is below 10^⁻²⁶, so the null hypothesis that all 43 algorithms perform equally is rejected overwhelmingly in all four cases and the observed ordering is not a chance result of sampling. The evidence is strongest on CEC-2017, the largest suite, where the statistic is more than three times that of the next suite.

Aggregating the pairwise tests give a consistent picture of PAO’s standing. On CEC-2017 PAO records 1150 wins, 56 ties and 12 losses (94.4%); on CEC-2022, 445/32/27 (88.3%); on CEC-2019, 362/31/27 (86.2%); and on Classical-23, 641/42/283 (66.4%). On the three modern CEC suites PAO holds a net-positive record against every competitor — all 42 on CEC-2017 and CEC-2019, and 41 of 42 on CEC-2022, where GOA is exactly level — and sweeps 16, 18 and 20 opponents respectively without conceding a single function. The Classical-23 figure is markedly lower, but as Section 5.1 established this reflects the tie compression at the global optimum that flattens that suite rather than a genuine weakness.

Figure 17 places PAO and the eleven next-best algorithms on a rank axis for each suite, with a bar of length CD anchored at PAO; any algorithm lying beyond that bar is significantly worse than PAO at α = 0.05. The diagrams confirm PAO’s leading position on all three CEC suites and its close second on Classical-23. They also make an important limitation visible. Because the comparison involves k = 43 algorithms, the post-hoc critical difference is wide — CD = 10.51 on CEC-2017, 12.00 on Classical-23, 16.61 on CEC-2022 and 18.20 on CEC-2019 — since it must correct for 42 simultaneous comparisons over as few as ten functions. A wide CD means that the post-hoc test is conservative: many mid-field algorithms cannot be separated from PAO at this significance level even though they rank far below it. Under the Holm procedure PAO is shown to be significantly better than 34 of 42 competitors on CEC-2017, 24 on CEC-2022, 19 on CEC-2019 and 13 on Classical-23. The correct reading is that the omnibus and pairwise Wilcoxon results carry the statistical weight, while the CD diagrams provide a conservative visual confirmation rather than the primary evidence.

5.8. Scalability Across Dimensions

A global optimizer is only useful if its behavior degrades gracefully as the search space grows. To assess this, PAO was run against itself at increasing dimensionality: on the thirteen scalable classical functions (F1–F13) at D ∈ {10, 20, 30, 50, 100}, and on the twelve CEC-2022 functions at the two dimensions for which that suite is defined, D ∈ {10, 20}. All other settings were held at the values of Section 4, so dimension is the only varying factor. Results are summarized in Table 14 and Figure 18, with the full per-function values in Appendix G.

At D = 10 PAO solves most of the classical suite to machine precision — nine of the twelve minimization functions reach a mean error below 10^⁻⁸, several at the 10^⁻²⁶ level. As the dimension rises this precision is progressively lost: six functions remain below 10^⁻⁸ at D = 20 and three at D = 30, while at D = 50 and above no function is solved to that tolerance. Crucially, however, the absolute errors at high dimension remain small in the context of the search domain: at D = 100 the mean error is 4.72 on the Sphere (F1), 5.00 on the Step function (F6), 0.35 on the noisy Quartic (F7) and 0.95 on Griewank (F11) — values that represent a good solution of a 100-dimensional problem rather than a failure. The functions that degrade most are the ill-conditioned Schwefel 1.2 (F3, rising to 5.4 × 10^³) and the Rosenbrock valley (F5, 1.8 × 10^³), whose difficulty is known to compound sharply with dimension.

The Schwefel 2.26 function (F8) must be read separately because its global optimum itself scales with dimension (f* = −418.98 D), so a raw comparison of means across dimensions is meaningless. Expressed as the percentage of the optimum attained, PAO reaches 99.7% at D = 10 and 98.3% at D = 20, but only 81.7% at D = 30, 63.8% at D = 50 and 49.7% at D = 100. This is the clearest instance of dimensional degradation in the study and is reported openly: on a deceptive landscape whose optimum recedes with dimension, PAO’s relative performance falls to about half the ideal at D = 100 (Table G3).

On the modern suite the picture is more favorable, because the composition and hybrid functions are hard at both dimensions. The median error ratio between D = 20 and D = 10 is only 2.6, and on three functions (F9, F10 and the near-solved F5) the error at D = 20 is equal to or lower than at D = 10. Doubling the dimension therefore costs PAO less than a single order of magnitude on a typical function of this suite — evidence that the astronomical schedule of Section 3 continues to allocate exploration effectively as space grows.

The computational cost of the added dimensionality is strikingly small. Averaged over the classical suite, one run takes 1.52 s at D = 10 and only 2.02 s at D = 100 — a factor of just 1.33 for a tenfold increase in dimension, far below the linear growth a dimension-dominated cost model would predict (Figure 14b). PAO’s runtime is therefore governed by its fixed per-iteration overhead rather than by the dimension of the problem, which is a favorable property for large-scale use and is examined further in Section 5.11.

5.9. Population-Size Sensitivity (N = 20, 30, 40, 50, 60)

The population size N is the principal control parameter of PAO, governing the balance between exploration and per-iteration cost. To fix its value on evidence rather than convention, PAO was run at five population sizes, N ∈ {10, 20, 30, 50, 60}, on the CEC-2019 “100-Digit Challenge” suite (ten functions), with every other setting held at its default. For each N the mean solution error and its standard deviation over the 30 runs were recorded on all ten functions, giving a joint picture of accuracy and robustness. To summarize across functions of very different scales, the five population sizes were ranked on each function (1 = lowest error) and the ranks averaged; the results are reported in Table 15 and Figure 19, with the complete per-function values in Appendix H.

Solution quality improves steadily as the population grows, and N = 60 is best. Its average rank across the ten functions is 1.55, ahead of N = 50 (2.15), N = 30 (2.75), N = 20 (3.75) and N = 10 (4.80) a strictly monotone ordering with no reversal. N = 60 delivers the lowest mean error on six of the ten functions, and on two more (F2, the inverse Hilbert matrix, and F3, the Lennard-Jones cluster) the result is independent of population size, so these are effectively tied rather than losses. Only F8 and F10 do a smaller size edge N = 60, and there by a margin well within the run-to-run noise. There is no sign of over-population degradation anywhere in the tested range.

The benefit of a larger population is concentrated on the hardest, most exploration-demanding functions. Moving from N = 10 to N = 60 reduces the mean error on the Chebyshev problem (F1) by roughly four to five orders of magnitude (from 1.09 × 10^⁹ to 3.14 × 10^⁴), on Rastrigin (F4) by a factor of about 23, and on the modified Schwefel function (F7) by a factor of three, whereas the already-easy functions (F5 Griewank, F6 Weierstrass) improve only slightly. This is the expected signature of population size: a larger swarm sustains diversity and postpones premature convergence, which is decisive on rugged and high-conditioning landscapes and largely irrelevant once a function is easy.

Critically, the accuracy gains do not come at the price of stability — they are accompanied by it. Ranked by standard deviation, N = 60 is again the best size (average robustness rank 1.95, lowest variance on six of the ten functions). The effect is most dramatic on F1, whose standard deviation collapses from 1.37 × 10^⁹ at N = 10 to 2.71 × 10^⁻¹¹ at N = 60 — PAO becomes deterministic on the Chebyshev problem at the larger population — while the variance on Griewank (F5) and Happy Cat (F9) falls to the order of 10^⁻³. Larger populations therefore make PAO both more accurate and more repeatable, exactly the combination wanted from a default setting.

The cost of a larger population is proportionally more number of function evaluations per iteration: N = 60 performs three times as many evaluations per iteration as N = 20 and six times as many as N = 10, so the wall-clock time scales linearly with N (absolute runtimes are reported in Section 5.11). For a general-purpose global optimizer, however, solution quality and reliability are the primary objectives, and the evidence here shows that N = 60 buys order-of-magnitude accuracy improvements and a collapse in variance on the hardest problems for only a linear increase in cost, while degrading no function measurably. Although diminishing returns are visible beyond N = 50 on the easy functions, N = 60 remains the best or tied-best size on every function and the best overall by a clear margin. Accordingly, N = 60 is adopted as the default population size throughout this paper (Section 5.1, Section 5.2, Section 5.3, Section 5.4, Section 5.5, Section 5.6, Section 5.7, Section 5.8 and Section 5.9): it maximizes both accuracy and robustness within the practical range, concentrates its advantage precisely on the hardest functions where population diversity is decisive, and carries no measurable penalty on any function in the suite.

5.10. Runtime: PAO’s Speed Rank and What it Buys

PAO is a computationally expensive optimizer, and this is stated clearly. Measured as mean wall-clock time per run and ranked against the other 42 algorithms, PAO places 40th of 43 on Classical-23, 39th on CEC-2017, 41st on CEC-2022 and 36th on CEC-2019 — a speed rank in the range 36–41, that is, among the slowest quarter of the field on every suite. In absolute terms one run costs 1.72 s, 1.92 s, 1.80 s and 2.25 s on the four suites respectively, against field medians of 0.28 s, 0.47 s, 0.33 s and 1.13 s, so PAO is roughly two to six times slower than the median competitor. These numbers are reported in Table 16 and Figure 20; the full per-algorithm timings are in Appendix I. No attempt is made to present this as anything other than a genuine cost of the method.

Three factors explain the ranking, and none of them is an implementation defect. First, PAO spends more function evaluations per iteration than its competitors. Every algorithm in the comparison was run with the same population size (N = 60) and the same 500 iterations, so the population is not the cause; the extra cost comes from the phase operators. On top of the base engine, which evaluates the whole population, an active prayer window evaluates a further subset, giving about 1.12 N evaluations per iteration (Eq. 22) instead of N — roughly 12% more work for the same iteration count. Second, each PAO iteration evaluates the astronomical position schedule of Section 3 — solar-altitude terms that determine the non-monotone exploration profile — which is arithmetically heavier than the single vector update used by simpler swarm methods. Third, the cost is dominated by this fixed per-iteration overhead rather than by problem dimension: as Section 5.9 showed, a tenfold increase in dimension raises the runtime only 1.33-fold. Normalized per function evaluation, PAO costs roughly 57–75 microseconds, and this figure is nearly independent of the problem size.

The essential question is not whether PAO is slow but whether the time is well spent, and on this the evidence is unambiguous. The fastest algorithms in the pool are also among the least accurate. HS, the fastest method on every suite (0.007–0.024 s per run, roughly 200 times quicker than PAO), ranks 42nd, 36th, 41st and 37th of 43 in accuracy on the four suites; DAOA, the second fastest, ranks last (43rd) on all four. PAO, by contrast, ranks first on CEC-2017, CEC-2019 and CEC-2022 and second on Classical-23. Figure 20b plots speed rank against accuracy rank on CEC-2017 and shows no systematic advantage to being fast: the cheap algorithms cluster in the low-accuracy region, and the additional time PAO spends is converted into solution quality rather than wasted. It is also worth noting that PAO is not the slowest method — several competitors are slower still, and on CEC-2019 the slowest (ZOA, 42.9 s) takes nineteen times longer than PAO while ranking well below it.

In absolute terms the cost is modest: under two and a half seconds for a complete 30 000 evaluation run on a standard desktop. For the great majority of engineering applications — including the constrained design problems of Section 6, where a single objective evaluation may involve a simulation far more expensive than the optimizer’s own arithmetic system. This overhead is negligible against the value of a better solution. The honest summary is therefore that PAO trades wall-clock speed for accuracy and robustness: it is among the slowest quarters of the field, for which it delivers the best or second-best solution quality on every suite tested. Where evaluation budgets are extremely tight or the objective function is itself trivially cheap, a faster method may be preferable; where solution quality governs, PAO’s cost is well justified.

5.11. Time Complexity

Where Section 5.10 reported wall-clock speed, this section addresses algorithmic complexity: how PAO’s cost grows with problem size and population, and how much work it performs relative to the other 42 algorithms. Two forms of evidence are given: the theoretical complexity derived from the algorithm’s structure, and the measured complexity from the standard CEC timing protocol. Results appear in Table 17 and Figure 21 and Figure 22, with full data in Appendix J.

For a run of T iterations with a population of N agents in D dimensions, and writing Cf for the cost of one objective evaluation, PAO requires

$$O( T \times N \times ( D + log N + Cf ) )$$

The three terms correspond to the position update over D coordinates for each agent, the ranking of the population that the prayer-time schedule requires at each iteration, and the objective evaluations. Every one of the 42 competitors falls in the class O(T × N × (D + Cf)), so PAO is the only algorithm in the pool carrying an additional log N term. The term is additive rather than multiplicative, so PAO stays in the same polynomial class as every standard swarm method — linear in the number of iterations, linear in the population and linear in the dimension. At the settings used throughout this paper (N = 60, log₂ N ≈ 5.91, D = 30) the sorting term raises the per-iteration arithmetic by roughly 19% over an equivalent method without it, and its relative weight falls as the dimension grows because D dominates log N (Figure 22b). PAO is therefore not for growth rate more expensive than its competitors in any meaningful sense; it carries a modest constant-factor and logarithmic overhead.

The empirical counterpart follows the standard convention. $T₀$ is the time for a fixed reference computation, measured here as 0.011449 s; $T₁$ is the time to evaluate the objective function alone for the full budget (mean 0.7786 s across the 23 classical functions); and $T₂$ is the mean time for the complete algorithm over five runs. The reported complexity is the dimensionless ratio $(T₂ - T₁) / T₀$, which isolates the optimizer’s own overhead from the cost of the objective function and normalizes it by machine speed, so the figure is comparable across hardware. Averaged over the 23 functions PAO records a measured complexity of 511.68, ranging from 415.8 to 613.4 with a median of 494.6 — a narrow band, showing that PAO’s overhead is stable and largely independent of the function, as the theory predicts.

On this measure PAO ranks 34th of the 40 algorithms for which a valid measurement exists (three — HS, DAOA and DOA — could not be ranked, as explained below). This differs from Section 5.11, where PAO’s wall-clock speed rank was 36–41 of 43, whereas its complexity rank is 34 of 40: speed measures total elapsed time including objective evaluations, whereas the complexity ratio strips those out. Both place PAO in the heavier part of the field, the complexity rank being the more favorable because much of PAO’s wall-clock cost comes from the evaluations implied by N = 60 rather than heavier per-agent arithmetic.

The ranking carries one qualification, which the data makes explicit. The complexity ratio is only meaningful when an algorithm consumes the full evaluation budget; if it terminates early, $T₂$falls below $T₁$ and the ratio becomes negative and uninterpretable. Such measurements were excluded. PAO was valid on all 23 functions, but several of the apparently lightest algorithms were not: CSA, ranked first at 1.55, is valid on only 7 of 23 functions, SFOA on 16, SHS on 17, MRSO on 18 and TSO on 19, while HS, DAOA and DOA produced no valid measurement on any function. Their low scores therefore partly reflect not using the budget rather than performing the same work more cheaply. Restricting the comparison to the 31 algorithms measured on all 23 functions, PAO ranks 25th of 31 — still in the heavier group, but the gap to the lightest fully-measured method (FOA, 17.31) is the honest comparison rather than the gap to CSA. Figure 22a plots measured complexity against the number of valid functions and made this pattern visible.

The complexity analysis supports the same conclusion as the runtime discussion, on firmer ground. PAO is not in a worse complexity class than its competitors: like every other method in the pool, it is linear in iterations, population and dimension, and its one structural addition — the ranking required by the prayer-time schedule — adds a log N term that is additive, modest at these population sizes, and less significant as the dimension grows. The measurements show a larger constant factor, consistent with the richer per-iteration update described in Section 3 and with the deliberate choice of N = 60 justified in Section 5.10. Because the cost is dominated by this fixed per-iteration overhead rather than by problem size, it grows slowly with dimensions: the 1.33-fold increase from D = 10 to D = 100 reported in Section 5.9 — which is the property that matters for scaling to larger problems. Set against the accuracy results of Section 5.2, Section 5.3 and Section 5.4, where PAO ranks first on three of the four suites, this overhead is a reasonable price and does not compromise the method’s asymptotic behavior.

5.12. Ablation Study

An ablation study isolates the contribution of each component of PAO by disabling it and re-running the complete evaluation. Nine configurations were compared on the CEC-2017 suite (29 functions, 30 independent runs, 500 iterations): the full algorithm; five variants each removing one of the prayer-time operators (Fajr, Eq. 16; Dhuhr, Eq. 17; Asr, Eq. 18; Maghrib, Eq. 19; Isha, Eq. 20); a BaseEngineOnly variant with all five operators and the schedule disabled; and two variants replacing the astronomical solar schedule with a fixed and a random schedule respectively. Every variant was assessed by per-function mean error, by the two-sided Wilcoxon rank-sum test against the full algorithm (α = 0.05), and by average rank across the 29 functions. Results are summarized in Table 18 and Figure 23 and Figure 24, with per-function values in Appendix K.

The five prayer-time operators show no measurable individual contribution. This is the clearest outcome of the study, and it is reported as measured. Removing any single operator leaves performance statistically unchanged on the great majority of functions. Removing Maghrib or Isha changes nothing significant on any of the 29 functions (0 wins, 29 ties, 0 losses in each case); removing Fajr likewise yields 0–29–0, and in fact reproduces the full algorithm’s mean exactly on 28 of the 29 functions, differing only on F20 and there by 0.09%. Removing Asr favors the full algorithm on just two functions and removing Dhuhr produces an even split (four functions each way). In median terms the penalty for deleting any one operator lies between 0.00% and 0.18%, which is indistinguishable from run-to-run noise. On the average rank metric several of the reduced variants — no-Fajr, no-Isha — in fact place above the full algorithm. The evidence therefore does not support attributing PAO’s benchmark performance to the individual prayer-time update rules.

The solar schedule is a variance-reduction mechanism, not a mean-improvement one. The schedule ablations give a more nuanced and, on inspection, more informative result. Judged on the typical function, replacing the astronomical schedule does not hurt and may help: the random-schedule variant is a median 9.7% better than the full algorithm, wins on 15 functions against 8 losses, and attains the best average rank of all nine configurations (4.00 against the full algorithm’s 4.97). The fixed-schedule variant is likewise a wash (11 wins, 13 losses) and ties the full algorithm exactly on average rank (4.9655 for both). Taken alone these figures would suggest the schedule is unnecessary. The mean behavior, however, tells the opposite story: the random schedule is on average 371% worse, the fixed schedule 315% worse and the base engine 208% worse. The reconciliation is that the ablated variants fail abruptly on a small number of functions. Randomizing the schedule degrades the mean error by a factor of 105 on F9 and by 523% on F6, while fixing it degrades F9 by 5807%, F6 by 1664% and F1 by 1537% — six functions in each case exceed a 50% penalty. The astronomical schedule therefore buys protection against worst-case failure on specific multimodal landscapes, at the cost of a small typical-case penalty. That is a defensible role for the component, but it is a different claim from the one usually made for such mechanisms, and it should be stated in those terms.

Read as a whole, the ablation supports a narrower claim than component-wise decomposition usually implies. It does not establish that each prayer-time operator earns its place: on this suite, at this population size and budget, four of the five are individually removable without a statistically detectable change, and one of them (Fajr) is very nearly inert. Nor does it establish that the astronomical schedule improves typical performance; on the median function a random schedule is at least as good. What it does establish is that the composite configuration is the most reliable one: every ablated variant, including those with better median behavior, incurs order-of-magnitude failures on a subset of functions that the full algorithm avoids. PAO’s benchmark standing in Section 5.2, Section 5.3 and Section 5.4 is built on consistently avoiding such failures, and in that specific sense the components are doing work — collectively and in the tail, rather than individually and on average.

6. Real-World Engineering and Application Problems

Benchmark suites establish an optimizer’s numerical behavior, but engineering practice adds constraints, mixed variable types and strongly non-convex feasible regions that synthetic functions do not reproduce. This section evaluates PAO on twenty-two constrained real-world design problems — among them springs, beams, trusses, a pressure vessel, gearboxes, bearings and a car side-impact structure — drawn from the standard metaheuristic test collections, against the same pool of 43 algorithms and under the same protocol used throughout Section 5. Background theory for each application — its engineering goal, design variables, objective function and constraints — together with a labelled diagram of every problem and the complete per-problem results, is given in Appendix L (Real-World Engineering Problems — Detailed Results). This section first sets out the constraint-handling method and the test suite, then compares the designs PAO obtains with the best values published in the literature, reports the same statistical tests used for the benchmark results, and identifies the problems on which PAO does not reach the published optimum.

6.1. Constraint-Handling Method

All twenty-two problems are inequality-constrained minimizations of the form $min f\left(x\right)$ subject to $gj\left(x\right) \le 0$. Constraints are enforced by the static penalty method of the SNS reference implementation [(Talatahari, Azizi, et al., 2021)], in which the fitness minimized by every algorithm is

$$F\left(x\right) = f\left(x\right) + \Sigma j \nu j · max(0, gj(x\left)\right)$$

where $\nu j$are the problem-specific violation factors supplied with each formulation. Only violated constraints are penalized, so a feasible design satisfies $F\left(x\right) = f\left(x\right)$ and the reported fitness of any feasible solution is directly comparable with published objective values. The penalty is identical for all 43 algorithms, so the comparison is not biased by constraint handling. Two implementation choices deserve mention. First, because these problems have non-uniform vector bounds — each design variable has its own limits — the search is conducted in the normalized cube [0, 1]^d and the wrapper maps $x = LB + u · (UB - LB)$ inside the objective. This allows published algorithm codes that assume scalar bounds to run unmodified, and since the same mapping is applied to every algorithm the comparison remains fair. Second, designs producing complex or non-finite objective values (possible through fractional powers of negative intermediates in several formulations) are assigned a large penalty, and integer or discrete variables (gear teeth, disc counts, standard sections) are rounded or mapped to their catalogue inside the objective so that any continuous optimizer can be applied.

6.2. Classical Constrained Design Problems

The first group comprises the classical mechanical and structural design problems that form the standard validation set in the metaheuristic literature (Kumar et al., 2020; Talatahari, Azizi, et al., 2021): the speed reducer (P01, 7 variables, 11 constraints), the tension/compression spring (P02), the pressure vessel (P03), the three-bar truss (P04), the gear-train ratio problem (P05, a pure integer formulation), the cantilever beam (P06), the I-beam deflection problem (P07), the tubular column (P08), the piston lever (P09), the car side-impact design (P11, 11 variables with discrete section choices), the welded beam (P12) and the reinforced-concrete beam (P13, mixed discrete–continuous). These span two to eleven design variables and up to eleven constraints, and their global optima are well established, which makes them a strict test of solution accuracy rather than of ranking alone. The corrugated bulkhead problem (P10) of the original collection is excluded from this study because its published formulation admits designs that are numerically degenerate under the shared penalty wrapper, which would make cross-algorithm comparison unreliable.

6.3. Domain Applications

The second group extends the evaluation to ten further applications that exercise different engineering domains and numerical characteristics (Himmelblau, 1972; Kumar et al., 2020): the multiple-disk clutch brake (P14) and rolling-element bearing (P15, a maximization of dynamic load capacity exposed as minimization of its negative, 10 variables and 9 constraints); the Belleville disc spring (P16) and hydrostatic thrust bearing (P17), both strongly non-linear; the robot gripper (P18), whose objective is a min–max force range evaluated over a sampled displacement trajectory; the gas transmission compressor (P19) and stepped cone pulley (P20, with three equality constraints encoded as |h| − tol); the planetary gear train (P21, all-integer teeth counts); the ten-bar planar truss (P22), where each evaluation requires a finite-element solve and 18 stress and displacement constraints are imposed; and the Himmelblau heat-exchanger design (P23). This group is deliberately harder and less standardized than the first: it introduces simulation-based objectives, equality constraints, integer programming and maximization, and therefore tests whether PAO’s benchmark behavior transfers to realistic engineering formulations.

6.4. Results and Discussion

Each of the 43 algorithms was run 30 independent times on every problem for 500 iterations with a population of N = 60 — the default justified on evidence in Section 5.10 — using a common random-seed sequence across algorithms. Averaged over twenty-two problems, PAO attains a mean rank of 2.955 and ranks first of the 43 algorithms, ahead of APO (4.455), GOA (6.841), SFOA (8.341) and PGA (8.886). The Friedman omnibus test rejects rank equality decisively (χ² = 725.38, df = 42, p ≈ 2.1 × 10^⁻¹²⁵). PAO ranks first or second on twelve of the twenty-two problems, within the top four on eighteen, and never falls below eighth. Aggregated over all 42 opponents, the Wilcoxon rank-sum test (α = 0.05) gives PAO 843 wins, 34 ties and 47 losses (91.2%), with a net-positive record against every one of the 42 competitors without exception, including APO (12–1–9), and it sweeps 19 opponents without conceding a single problem.

For engineering problems, the decisive question is not ranking but whether the optimizer recovers the best design known in literature. PAO matches or improves the reference optimum on eighteen of the twenty-two problems, typically to eight or more significant figures — for example 2994.4245 on the speed reducer, 0.0126652 on the spring, 6059.7143 on the pressure vessel, 1.7248523 on the welded beam and 263.8958 on the three-bar truss (Table 19). On three problems it returns a better objective than the listed in the table reference: the rolling-element bearing (−82 842 against −81 859), the gas-transmission compressor (2.293 × 10^⁶ against 2.965 × 10^⁶) and the stepped cone pulley (16.088 against 16.63). These three belong to the group whose reference values the source implementation itself describes as approximate targets rather than exact optima, so they are reported as improvements on the stated reference under this formulation rather than as new global optima. Considered across the whole pool, PAO obtains the best design found by any of the 43 algorithms on twenty of the twenty-two problems.

Four problems are not solved to their tabulated optimum: the gear train (P05), the clutch brake (P14), the hydrostatic thrust bearing (P17) and the robot gripper (P18). This is reported plainly, together with the diagnostic that explains it: on all four, none of the 43 algorithms reaches the tabulated value, and on three of them PAO is the best or joint-best method in the pool. On P17 PAO attains 19 505 where the field best is also 19 505 and the reference is 1 625; on P14 it is joint-best at 0.3137 against a reference of 0.2352; on P05 it returns 2.70086 × 10^⁻¹², which is the classical literature optimum, while the reference value stored in the source code carries a −10^⁻⁴ offset that no algorithm can attain. Only on the robot gripper (P18) does another method finish ahead of PAO (HBA, 4.288 against 4.324), and there too the reference of 2.529 is unreached by every algorithm. The pattern indicates discrepancies between the coded formulations and the published reference values — a possibility the source implementation explicitly flags for problems 15, 16, 18, 20 and 21 — rather than a failure of the optimizer. PAO’s weakest ranking on the suite is the piston lever (P09, eighth of 43) — the only problem on which a group of competitors finds better mean designs, and even there PAO recovers the reference optimum 8.41270 in its best run.

PAO’s advantage is in mean performance across runs and not merely in its best run: its per-problem ranks are computed on the 30-run mean, so a first or second place indicates that the algorithm reaches near-optimal designs reliably rather than occasionally. The computational cost follows the pattern established in Section 5.11 and Section 5.12: PAO averages 1.69 s per run against a field median of 0.67 s, a speed rank of 40 of 43. As argued there, for engineering design this overhead is immaterial — a single finite-element evaluation in the ten-bar truss problem already dominates the optimizer’s own arithmetic, and two seconds of search time is negligible against the value of a better structural design.

On twenty-two constrained real-world design problems spanning mechanical, structural, civil, automotive, robotic and process engineering, PAO ranks first of 43 algorithms with a mean rank of 2.955, recovers or improves the reference optimum on eighteen, obtains the best design in the entire pool on twenty, and holds a net-positive statistical record against all forty-two of its competitors. Taken with the benchmark results of Section 5.2, Section 5.3 and Section 5.4, this establishes that PAO’s numerical performance transfers to constrained engineering practice, which is the setting these algorithms are intended for, (see Table 20).

Figure 25. PAO on the twenty-two constrained engineering design problems. (a) Average rank of the twelve leading algorithms across the problems (lower is better); PAO ranks first. (b) PAO’s rank on each individual problem, with its average rank marked; PAO is first or second on twelve problems, within the top four on eighteen, and never falls below eighth (piston lever, P09).

Figure 25. PAO on the twenty-two constrained engineering design problems. (a) Average rank of the twelve leading algorithms across the problems (lower is better); PAO ranks first. (b) PAO’s rank on each individual problem, with its average rank marked; PAO is first or second on twelve problems, within the top four on eighteen, and never falls below eighth (piston lever, P09).

7. Discussion and Limitations

This section steps back and reviews the evidence honestly. It identifies cases where PAO does not come first and explains the reasons, including test problems that are too easy to separate the algorithms and problems where the constraint-handling method is the limiting factor. It also discusses the extra computing time PAO needs compared with simpler methods, how sensitive the results are to the chosen settings, and which conclusions could be affected by the way the experiments were designed. The purpose is to state clearly what has been shown, what has not, and what should be done next.

7.1. Where PAO does not Win, and why

PAO leads four of the five test collections, and an honest account of the fifth and of the individual functions it loses is more informative than the aggregate. Three distinct situations arise.

(i)

Rank compression on easy suits. On the classical 23-function set PAO places second (mean rank 13.152) behind SOA (11.261). Section 5.1 traced this to the structure of the metric rather than to search quality: on the easy separable functions a large fraction of the 43 algorithms drive the mean exactly onto the optimum, so under fractional tie-breaking a block of k tied algorithms all receive the single rank (k+1)/2 — up to eighteen algorithms share one rank on F11. PAO is penalized twice by this: on F1, F9, F10 and F11 its mean sits marginally outside the tie block and it is displaced below the entire block, and on F16, F17 and F19 it lies inside the block and its accuracy is not distinguished from that of thirteen to sixteen peers. The interpretation is supported by what happens when the compression disappears: on the shifted and rotated CEC suites, where ties at the optimum vanish, PAO ranks first while SOA falls to 32nd and CSA to 42nd on CEC-2017. The classical result is therefore a property of an insufficiently discriminating benchmark, not evidence of weaker search.

(ii)

A genuine near-tie at the top. On CEC-2022 PAO’s first place is narrow. Its mean rank of 4.375 is only marginally better than GOA’s 4.667, and the head-to-head Wilcoxon record is exactly level (5 wins, 2 ties, 5 losses). The two methods are best described as sharing the leading position on that suite rather than separated by it, and the paper states this rather than claiming a clean win.

(iii)

Specific landscape weaknesses. PAO has a small number of reproducible failure modes, concentrated on deceptive and fractal landscapes. The clearest is the Schwefel function, whose global basin lies far from the strong local basins: PAO ranks 14th of 43 on CEC-2017 F10, and under increasing dimension its attainment on the classical Schwefel function falls from 99.7% of the optimum at D = 10 to 49.7% at D = 100. The second is the Weierstrass function (CEC-2019 F6, rank 11), whose fractal roughness frustrates late-stage refinement. Others are the Rosenbrock valley on CEC-2022 (F2, rank 11) and two of its composition functions, and the piston lever among the engineering problems (P09, eighth of 43). In each case the difficulty is a landscape feature that specifically defeats the late refinement phase of the astronomical schedule, and each is narrow enough to be stated precisely rather than absorbed into an average.

7.2. Runtime and Computational Overhead

PAO is an expensive optimizer, and this is not disputed. Its wall-clock speed rank is 36–41 of 43 across the four benchmark suites and 40 of 43 on the engineering problems, with one run costing 1.7–2.3 s against field medians of 0.3–1.1 s. Three factors explain the cost. The population size N = 60 selected in Section 5.10 performs twice the objective evaluations of a method using N = 30; each iteration evaluates the solar-altitude schedule, which is arithmetically heavier than a single vector update; and the population must be ranked each iteration, contributing the log N term in the complexity expression O(T·N·(D + log N + Cf)) — the only such term in the pool, and an additive one that leaves PAO in the same polynomial class as every competitor. Two mitigating observations follow from the measurements. First, the cost is dominated by fixed per-iteration overhead rather than by problem size: a tenfold increase in dimension raises runtime only 1.33-fold, so the overhead does not compound as problems grow. Second, cheapness in this pool does not buy quality — HS, the fastest method, is roughly 200 times quicker than PAO yet ranks 36th to 42nd in accuracy, and DAOA ranks last on all four suites. The overhead is therefore a deliberate exchange of time for solution quality, and it is immaterial wherever a single objective evaluation involves a simulation, as in the finite-element truss problem of Section 6.

7.3. Parameter Sensitivity

PAO exposes one principal control parameter, the population size. The study of Section 5.10 across N ∈ {10, 20, 30, 50, 60} on CEC-2019 found a strictly monotone improvement in averaged rank (4.80 → 3.75 → 2.75 → 2.15 → 1.55) with no reversal and no sign of over-population degradation, accompanied by a parallel improvement in robustness. The gains concentrate on the hardest functions — four to five orders of magnitude on the Chebyshev problem — and are negligible on easy ones, which is the expected signature of population diversity. This behavior is convenient in practice, because it means performance is not sensitive to fine tuning: any sufficiently large population performs well, and N = 60 was adopted because it was best on every criterion tested within a practical range. The limitation is that the range itself was bounded at 60 and examined on a single suite; the point of diminishing returns beyond that value, and any interaction between N and the iteration budget, remain uncharacterized. The remaining schedule constants are fixed by the astronomical model rather than tuned, which removes them as free parameters but also means they were not optimized.

7.4. Threats to Validity

Section 5.13 does not support a component-wise decomposition of PAO’s performance. Removing any one of the five prayer-time operators leaves results statistically unchanged on the great majority of the 29 CEC-2017 functions — Maghrib, Isha and Fajr each produce 0 wins, 29 ties and 0 losses — and the no-Fajr variant reproduces the full algorithm’s mean exactly on 28 of 29 functions. On the average rank metric a random schedule ranks ahead of the full configuration. What the ablation does support is narrower and should be claimed instead: the ablated variants fail catastrophically on roughly six functions each, by factors of up to 10^², which the full configuration avoids, so the astronomical schedule behaves as a variance-reduction and tail-protection mechanism rather than a mean-improvement one. Two consequences follow. The claim that each operator earns its place is not evidenced by these experiments; and the near-identity of the no-Fajr variant is strong enough to warrant verifying that the Fajr branch is exercised as intended under the reported settings. The ablation also covers a single suite at one dimension and tests only single-operator removals, so redundancy in isolation does not exclude necessity in combination.

Comparing 43 algorithms over as few as ten functions makes the post-hoc critical difference wide (10.5–18.2), so many mid-field algorithms cannot be separated from PAO at α = 0.05 even when they rank far below it. The Friedman and pairwise Wilcoxon results carry the statistical weight here; the CD diagrams are a conservative visual confirmation and should not be read as evidence of equivalence. All conclusions are also conditional on the shared experimental budget — 30 runs, 500 iterations, N = 60 — and a different budget could reorder closely spaced competitors.

All algorithms were run from published implementations on one machine under a single MATLAB version, so results inherit whatever implementation quality those codes carry; a competitor disadvantaged by its own published code would be disadvantaged here. In the engineering study, four problems (P05, P14, P17, P18) are not solved to their tabulated optimum by any of the 43 algorithms, which indicates discrepancies between the coded formulations and the cited reference values rather than optimizer failure; the source implementation itself flags several of these values as approximate. Cross-algorithm comparison remains valid because every algorithm optimizes the identical coded formulation under the same penalty wrapper, but absolute claims against the literature for those four problems should be treated cautiously. Constraints were handled by a single static-penalty scheme, and a different constraint-handling technique could alter the engineering ranking.

Finally, the no-free-lunch theorem (Wolpert & Macready, 1997) guarantees that no optimizer is superior across all problem classes. The evidence here establishes that PAO performs strongly on the shifted, rotated, hybrid and composition landscapes of the modern CEC suites and on constrained engineering design; it does not establish superiority in general, and Section 7.1 identifies the landscape families on which it is beaten.

8. Conclusion and Future Work

This paper introduced Pray Algorithm Optimization, a population-based metaheuristic whose exploration–exploitation schedule is derived from solar geometry rather than from a monotone decay law, producing a non-monotone profile in which diversity is restored at fixed, planned points of the run. PAO was evaluated against 42 competing algorithms over 75 benchmark functions and 22 constrained engineering design problems, with 30 independent runs of 500 iterations throughout. The numerical outcome is summarized in Table 21. PAO ranks first of 43 algorithms on CEC-2017, CEC-2019, CEC-2022 and the engineering suite, and second on the classical set, where tie compression at the global optimum limits the resolution of the metric. Friedman tests reject rank equality on every collection with p < 10^⁻²⁶, and PAO holds a net-positive Wilcoxon record against 41 or 42 of its 42 competitors on each of the modern suites. On the engineering problems it recovers or improves the reference optimum on 18 of 22 designs and obtains the best solution found by any algorithm on 20. Diversity measurements show the balance is adaptive, ranging from 4% exploration on unimodal landscapes to 52% on the deceptive Schwefel function, and the cost of this behavior is a speed rank of 36–41 of 43 with no increase in asymptotic complexity class.

PAO is defined as continuous domains. A binary formulation — mapping the continuous position through a transfer function, as is standard for feature selection, knapsack and scheduling problems (Mirjalili & Lewis, 2013) — would extend it to combinatorial tasks, and the non-monotone schedule is a plausible advantage there because binary landscapes are highly multimodal. A multi-objective version using Pareto dominance with an external archive and crowding-distance selection (Deb et al., 2002) is a natural companion, since most engineering design problems of Section 6 have competing objectives (weight against cost, capacity against reliability) that were collapsed into a single measure here.

The failure modes identified in Section 7.1 suggest specific hybridizations rather than generic ones. Because PAO’s weaknesses lie in late-stage refinement on deceptive and fractal landscapes, coupling it with a local search — a memetic scheme applying pattern search or a quasi-Newton step to the incumbent — targets exactly the deficit observed on Schwefel and Weierstrass. Given the ablation result of Section 5.13, an equally important line is mechanism refinement: verifying and, if necessary, redesigning the individual operators so that each demonstrably contributes, and running the intermediate configurations (operator pairs and triples) needed to separate redundancy in isolation from necessity in combination.

PAO’s cost is dominated by fixed per-iteration work and scales almost flat in dimension, which makes it a good candidate for parallelization. The operators act on agents independently within an iteration, so objective evaluations map naturally onto multiple cores or a GPU, and an island model with periodic migration would suit the larger populations that Section 5.10 shows to be beneficial. This is the most direct route to remove the runtime disadvantage documented in Section 7.2.

Beyond the design problems studied here, the natural targets are those where objective evaluation is expensive and gradients are unavailable — neural architecture and hyperparameter search, feature selection in high-dimensional data, energy-system scheduling, and structural optimization at larger scale. Two methodological extensions would strengthen the approach itself: an adaptive population size that grows or contracts on measured diversity rather than being fixed in advance, and a constraint-handling comparison against feasibility rules and ε-constrained methods to establish how much of the engineering performance depends on the static-penalty choice. Evaluation at higher dimensions than those reported here would also test whether the graceful degradation observed up to D = 100 continues at larger scale.